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ABSTRACT 


A study is made of the problem of periodic heat flow to the 
ground in a freeze-thaw situation. Analytic and numerical models are 
developed for the purpose of studying the effects of various parameters 
on the ground thermal regime. In particular, the effects of surface 
cover, soil properties and meteorological parameters on the average 
ground temperature below the active layer are studied. The main emphasis 
has been toward obtaining an indication of the effects of the parameters 
rather than an overall predictive model. 

The results obtained indicate the relative importance of the 
various parameters. They also indicate that under certain conditions 
simplifying approximations can be made to the environmental boundary 
conditions. Analytical methods of estimating the quantitative effects 
of surface disturbances are proposed. 


The results are in general agreement with field observations. 
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PART I 
INTRODUCTION 


The problem of heat flow into a thawing and freezing soil has 
been studied rather extensively in the past especially by those countries 
close to the Arctic Region namely, Russia, Canada, United States and 
the Scandinavian countries. Despite these efforts, an exact mathe- 
matical solution to any practical thermal soil problem is still 
unavailable. General solution of the temperature distribution in the 
soil involves analysis of the transcient heat flow in the soil which is 
affected by the thermal properties, namely, the specific heat, latent 
heat and conductivity, all of which are temperature-dependent. There- 
fore in order to find an exact mathematical solution, the specific 
problem must first be idealized. Problems of this nature are boundary 
value problemS and therefore accurate specific boundary conditions must 


be applied to the surface where the heat flow is known. 
1.1 Review of some analytical techniques 


The first mathematical solution obtained for this problem is 
by Neumann (1) but was later adapted to the calculation of frost depth 
penetration in soil by Berggren (2). The first published literature on 
the rate of ice formation in still water was by Stefan (3) and thus 
hence originate the title "Stefan Problem". Most of the theoretical 
frost depth equations used even today can be shown to be a modified form 
or a variation of the "Berggren formula" (or Neumann) even the "Stefan 
formula". Later Aldrich and Paynter (4, 5) presented the "Berggren 


formula" in a modified form so that it could be adapted to layered 
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"modified Berggren formula" has been adopted 


systems. Since then the 
by the Corps of Engineers and the Highway Research Board in the United 
States with some modifications to accommodate the specific conditions. 

Consider now Neumann's solution applied to a semi infinite 
soil initially at constant temperature, TQ: Then assume that the 
surface temperature drops below freezing and remains so for all t > 0. 
There exists now two phases, the frozen and the unfrozen region. 


The temperature distributions for the frozen and the unfrozen 


regions are given by: 


ey =e Cle/eren) erf (y/2(a,t)*/*) (1.1) 
Te = T, ~{1@_-2,)/erfe Maga)! ?] erfc(y/2(a,,.t)"/*)} 
(1.2) 


where Ty is the freezing temperature of the liquid and 


r is the proportional constant which is determined from: 


: 1/ 


2 
(emi /ee eae {Uk ¢ (og) ~(ag/aus) > 1/ 


2 
(T.-T,) e 


1/2 1/2 


= (ALYn)e,T ile 
[k, (a, ¢) Ty erfc Map/a¢) iB; OL 1) Ce 2 ( ) 
The interface position y is obtained from 


ve ar(a,t)/? aby 


Apart from Neumann (1, 6) and Stefan (3), the only two known 
exact solutions available, there are others solutions but they do not 


give the exact closed form solutions. The Variational Technique was 
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used by Biot and Daughaday (18) to study the heat conduction in a 
melting semi infinite solid with constant properties. Goodman (19) 
used the Heat-Balance Integral Technique to solve cases of melting 

in a semi infinite solid with a fixed boundary temperature and with a 
given heat flux at the boundary. The Heat-Balance Integral Technique 
changes the energy equation from a partial differential equation into 


an ordinary differential equation. 
x oT x 
i oc 5 OX = / Vik. Vivrdx Ci) 
oO oO 


The temperature profile is then assumed to be of a general polynomial 
form which is substituted into the integral equation. Integration 
produces an ordinary differential equation with time t as an 
independent variable. Poot (21) used the integral method to study a 
moving two-dimensional solidification problem. The complexity of this 
heat transfer problem is increased if the latent heat effect no longer 
occurs at a certain single temperature, but over a temperature range. 
The problem was described by Rubinshtein (22) and Weiner (23), both 
considering the soil undergoing "'n'" phase changes. In 1967 Tien and 
Geger (27) obtained a solution for the above problem using the Neumann 


approach by modifying the boundary conditions. 
1.2 Review of some of the numerical techniques used 


As reviewed previously, there are relatively few analytical 
solutions available for the practical heat transfer problems involving 
a change in phase of the soil due to freezing and thawing. In order 


to obtain a more general solution, other techniques have to be adopted. 
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Numerical methods are the most powerful and accurate of all the 
techniques available. 

Many of the existing numerical solutions to the practical 
heat transfer problem involving change of phase use a continuous 
description on the movement of the fusion front as well as a description 
of the temperature distributions. Modifying the conventional finite 
difference procedure, various authors developed different ways of 
solving the problem. Cochran (24) developed a simplified analysis 
by using a single lump of variable thickness to represent the frozen 
region. While the technique works the inaccuracies resulted from the 
lumping severely hindered its use. Otis (25) extended the works of 
Lightfoot by using the concept of a moving heat source to account for 
the latent heat. Forster (26) modified the conventional finite 
difference equation for fixed grid so that the travel of the fusion 
front) ss COneAmucdeky computed, but the method is limited for soil 
anitially or finally at the fusion temperature only. Murray and 
Landis (12) later developed two numerical approaches with a wider range 
of applicability and a greater degree of accuracy. One approach 
utilizes the variable space network where the number of grid points 
in each phase (frozen and unfrozen) ace constant. The second method 
was a refinement of the Forster work, it lead to a more accurate 
prediction of the temperature distribution and it also continuously 
computed the fusion front movement. 

Most of the above mentioned methods have difficulty in 


handling the situation where a fusion front overtakes a thawing front 
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or visa versa. This occurs in the case where the surface temperature 
is periodic and its range includes 0°C. The numerical model presented 


in this thesis will be free of this limitation. 
1.3 Recent interest 


Building of engineering structures such as oil and gas pipelines, 
roads and etc. in permafrost regions such as Alaska and Northern Canada, 
has been and still is a big problem for any contractor especially now 
with the present development and discovery of large energy resources. 
With the rapid growth of development, a growing concern for the arctic 
environment has caused both the government and industry to increase their 
rate of research on developing suitable methods for constructing and 
operating in permafrost areas. A considerable amount of work has been, 
and continues to be done in the area of developing comprehensive models 
of heat flux in freezing and thawing soils. Numerical methods have 
been, and continue to be developed by universities, government and 
industry to model ground thermal regime. The choice of numerical methods 
used depends on the problem under consideration. For one-dimensional 
analysis, finite difference nennade are in common use. The finite 
difference methods used are: ordinary forward difference, backward 
difference and the alternating direction. The forward difference has 
been used by many authors whereas the backward difference has been 
rarely used. The alternating-direction explicit methods were employed 
by Lachenbruch (13) and Doherty (14) while Fleming (15) used the 
alternating-direction implicit methods for solving the problem. The 


recent trend is towards the use of finite element methods, as were used 
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by Hwang ital (16) and Charlwood (17). These have been designed for 


engineering predictions of temperatures around engineering structures. 
1.4 Proposed program 


The purpose of this thesis will be to assess the effects of 
each of a number of parameters in this very complicated problem of 
periodic heat transfer to the ground in a freeze and thaw situation. 
This will be done by breaking the problem into three subproblems. They 
are: 

de Ground temperature distributions with imposed surface 
temperature. Here, the main objective will be to look at 
the effects that the various individual soil parameters 
have on the ground temperatures and the active layer. 

II. Ground temperature distributions with an insulating layer 
and an imposed insulating layer surface temperature. Here 
the main objective will be to look at the effects the 
insulating layer have on the ground temperatures and the 
active layer. 

III. As in (II) but with the environmental heat fluxes at an 
insulating layer surface boundary. Here, the main objective 
will be to look at the effects the main parameters have on 


the various temperatures. 
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PART LE 
IMPOSED SURFACE TEMPERATURE 


2.1 Analytical model for the ground temperature distributions with 


imposed surface temperature 


Die ed: Formulation of the model 


Assume the ground to be a semi infinite domain in which two 
phases exist as shown in Figure 1. Soil properties will be assumed 
to be uniform throughout a given phase, however conductivity, specific 


heat and the temperature gradient will be different for the two phases. 


Unfrozen 
Sousli 


Frozen 
Soil 


| Depth 


FIGURE 1 
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8. 
Applying the classical theory of heat transfer to the above conditions, 
the temperature distribution within a given phase may be calculated. 
Using the Fourier approximation for one-dimensional heat 


conduction and the first law of thermodynamics, the equation for the 


temperature Tut in the unfrozen region becomes 
0 (977 oes = OT fot (25) 
uf uf uf a 
where Oe Ke ong: 


At the moving interface between the frozen and the unfrozen 
regions, the heat flow into the interface’ equal to the heat flow 
out plus the latent heat of fusion associated with the interface motion. 


The conditions will be satisfied by the equations 


K ig OT g/dx + pL(ay/at) = k 


uf oT / ax ) 


be 


where y is the position of the interface. 


In the frozen region the diffusion equation is 
a .(8°T,/3x7) = 3T,/at (2.3) 
£ f £ 


The boundary condition at the ground surface will be assumed to be that 


of an imposed temperature varying sinusoidally during the year, that is 
T= T+ Aiesin re 
s Ss 


where T, is the average temperature and 


AT is the amplitude of variation 
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The boundary condition at x > ~ is 0T/dx = 0, that is zero heat flux. 


Non-dimensionalizing equations (2.1), 


equations 


and 


where 6 


become 
2 2 
r) 8 6/98 = Sp 28 OD) 
(k g/ke) 20,,/80 = - an/ot + 36 -/ 3% 


376 _/at- + $,(06,/23¢) 


= (T-T,)/AT 
GC = xX 
n = Vie. 
t =t/tec (anit, = ] year) 
Se = pc,AT/pL 
2 
aside Net, 

4 - 1/2 
where AS (ket AT/pL) 
similarly 
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The boundary conditions in non-dimensional form becomes 


at t= 0 § = 8. sp (Slatin Ware 
and 


at 0 > 90/9 = 0 


For the purposes of obtaining an approximate analytical solution, the 
above problem can be simplified down to "Stefan Problem" which is the 
simplest form of approximation. It will be used to predict the frost 
penetration depth and the average ground temperature. In this 
simplification it is assumed that the volumetric heat capacities of 

both the frozen and the unfrozen soil are negligible with respect to 
the latent heat of fusion of the soil water. This condition is sipeered 
when the Stefan Number (cAT/L) is small. When this condition is 


satisfied, the three governing equations become respectively 


2 2 
- = he 
P) 8 g/d 0 (247) 
= = + ) r) Dao 
(k ¢/Ke) 36. /8s dn/oat + 9 ¢/ @ (2.8) 
a6. /az = 0 | (2.9) 


Equations (2.7) and (2.9) with the unfrozen and the frozen phases can 


thus be integrated to give 


= (2.10) 
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that is within a given phase the temperature profiles are linear. 
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FIGURE 2 
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For the assumedform of the surface temperature, three different time 


domains may exist throughout the year (Figures 2 and 3). 
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(iii) 
FIGURE 3 


on 
> 
op 
- 
+ ‘© 
4 
£ 
5 
= 


fire’ winds . rte 
ey “hie ins ee 
‘ ns 


4 ‘ 


Mm 
Foi 
LY 
4 
4 > 
re 
. 
= 
y' 
: 
3 
4 
i 
3 
) 
& a 
' 
v 
} 
ri 
x 
. 
} 


A) 


a 
a 
t 


i 


7 
| 


¥ 


si 


Applying the boundary conditions to equation (2.10) and (2.11) 


Considering first equation (2.10) 


Of = C, a3 C., iG 
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; oy ‘ Ze 
at C= tN, ore = Nae 
C, = 2X 6./n 
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hee 6 (1-c/n) | (2.12) 


Considering now equation (2.11) 
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implying Cre =i 0 


hence Se a 8) (27215)) 
Substituting equations (2.12) and (2.13) into equation (2.8), it becomes 
- (kK p/kp) 0 ./n = - dn/dt (2.14) 


Equation (2.14) is the governing equation for this simplified problem. 
Drliene Solution 


Considering first region (1) and Figure 2(i), the boundary 


conditions are 


at n= 0, T= 


and» at n 
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The governing equation is 
(k g/Kp) 8. dt = ndn 
but -@ =6 + sin 2mT 
s s 
where os = (I -T,)/AT 
Integrating equation (2.14), it gives 


le = (k g/ Ke) {-cos 2nt/2a + ot} + const 
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from the boundary conditions, when n = 0, T = 
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From Figure 2 (i) 
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=~ (1-8,°) 


Substituting the above expressions into equation (2.15), it becomes 
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Equation (2.16) gives the depth of the thaw penetration. Considering 


now region (2) of Figure 2. 


ed 
rt 


wt 


(iit Bor oer paar + ACL ty ee at connie ae beh chi 
Ve 1 ; 3 ty ai ; ; 


gripe: « wy # . “Wa 


ioe 
j 
A wn 
1.7 g eth: 
a 4 yf 
1 
( » Ti © Air 
a 
*+ 
"7 ‘ 
z hanes 
f i rp. § 4 , ie 
i | ae - i 
& 
i 3 
7 
q - 
ra) ¥ 
| 
> 5 
i 
~ 
~ * 
Sy 
= 
, Ad 
‘ 
(Uk = 
; iY ye. 
Shak! hab NM: yee © , 
a TA Vers - = ¥ 


soya : , ‘ vi 
a ‘ x VFA 


OMNES GQ iNde > ataaihh oot wyolen dais avila gill aa eee aia? 
, : 7 us % Hi - ; 
; ‘ 7A 


a _ ss ‘ 


ts. ] 7 : é a , % j > ae 7 fh i — ; 
Paes veka Ghai oeAy ats Luo Deg piel LORRY eM a 
' ' h ny 7 4 § ef] = ( a ‘ atte ™? i . 1 
te WEA AG (2) ORG kg ay (Mae 7 
i iz ; - ¥ J ae 7 is 

el out id, ee ac 

| ‘ 5 yy 
- eT, 


16. 


aS c 
oe Ne /2 rien J ps dt 
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from region (1), 
2 “b 
i ee Hon he f @ dt (2.18) 
be 
7h Ss Pea Me 
but n /2 = Ne /2 
Te "b 
JESSY We ee cle Gna (2.19) 
iE is 
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Considering all the 3 regions 
a Tat] 
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Substituting equation (2.21) into equation (2.19), it gives 


— ae re — 

Gudeechje (Oyhdr suf poolsatan6 

gZ s ; Ss Ss 
"a b 


1 a \ 4 : ; ; a) ae 


me . : 
= > q /< 
‘ 4 #E 
‘ > Gis ts A ane Tt lhe ( 
ri Lo es) ot 4 > 
f 
' 
~ Aoy > = ’ 
— 


SA0t ah » (MESS) be pe CIS (i, 


ee on. ie eee 


. ~ ew 
wy ae % ¢ i. = ae aye S417 vas 
: if ; 
% fi ~ 
= 
+ ) - Noel 
a ae P , 7) 
j ; ' 
ioe 
° 
Th 
® 
i . 
! U r; 
ta = 
$ +o 
7 ye 
~» a9 
be . ne i 
a -_— a iF A 
= : i sh - 
‘oa * a a = Mire uJ 


a7 s 


a Peorifge bt ihinCandtecisl dgeel dt 
- a 
ON 
Serilak (ke) [omen da 
AG 
a 
but from region (1), 
"Dp — 2.1/2 a -1,= 
Sele sordtg= = (1R8 eye 6 (1/2-Binis (-6. ¥/r) 
a { i 
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+ (87)? pn} (2.22) 


Equation (2.22) gives the average ground temperature for the 


case when T < T.. 
s F 


Solutions to the above problem are for the case when the 


average temperature ine is less than the freezing temperature T,,. Using 
=) 


= 


the same approach as above, the solutions for the case when the average 


™m 


temperature Ty is greater than the freezing temperature Ty are 


: 241-8 7) + 0 [1/2-sin (8, )/n]} (2.23) 
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= > Le ie ies = Oe 
(0,-0,) = Arke/kg){ O[1/2-sin “ (-8,)/n + (1-8,°)"'/)} 


(2.24) 


Equation (2.23) and (2.24) gives the frost penetration and the 


average ground temperature for the case when T, > Ty These analytic 


results will be compared with the numericai results later in this part 


of the thesis. 
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yd ee Limitations of this analytic model 


» The assumptions that the properties of the soil are uniform 
is not true in many cases. This is because of the different 
soil textures, moisture content, soil compositions and other 
soil properties at different regions in the soil. All these 
properties contribute to the volumetric heat capacities, 
latent heat of fusion and the thermal conductivities of the 
soil. 

ii. Though the air temperature is not a true periodic function of 
time owning to the variability of the weather, it has many 
features in common with a periodic function, such as the 
fluctuations caused by the succession of day and night or 


winter and summer. 
2.2 Numericat Model 
Phe Dee Formulation of the model 


The main difference between the proposed method and most of 
the existing methods will be, instead of continuously computing the travel 
of the fusion front, this proposed method will keep track of the enthalpy 


h, as well as the temperature T, at any time t. Governing Equations: 


gh/at = k(9-T/dx-) (2.25) 

where h/pL = (T-T,)ec,/eL T < Ty 
(2.26) 
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The boundary conditions are 


C27) 
ye dT/dx = 0 


At T= Tp» the enthalpy is not a single valued function of temperatures, 
thus equation (2.26) cannot be substituted directly into equation (2.25) 
as in the case where no phase change occurs. This difficulty can be 
overcome either by keeping track of the position of the fusion fronts 
or by keeping track of the enthalpy. The later method will be used. 

Equation (2.26) is true when the entire phase change occurs 
at a single temperature. If the phase change happens to occur over a 
temperature range, equations (2.25) and (2.26) will be modified to take 
into account the unfrozen water content. The modifications will be 


presented at the end of this section. 
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Equation (2.26) becomes 


> = (pc ,AT/pL)6 O< 20 
o = (pc AT/pL)O + 1 8 > 0 
ie: > = 5,6 Gr<n0 
(2.23) 
>? = See aPadl : oy S10 
where Se = (ec ,AT)/pL 
and Sif = (pc, AT) /oL 
Equation (2.25) becomes 
(pL/t) d¢/9T = (KAT/x,”) se/ace 
on et 2 Z 2 
- do/aTt = [kt AT/pLx, ] 9 6/9z (2.29) 
For the frozen and unfrozen regions, the respective equations are: 
2 Z 2 
do/ot = [kt AT/pLx, ] 9° 8/9z (2.30) 
and d¢/ot = [k .t AT/pLx ay se /are (2.31) 
uf “c € 


choosing 


a5. 
[ket AT/pLx. }=1 
then equations (2.30) and (2.31) become 


d¢/dTt = 37e/ar- (2532) 
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2. 2 
and do/at = (Kg /ke) 0 6/93t (2.33) 
2 q 7 2 
ie: do/at = k'(90 6/9t°) Qs3h) 
where k= 6 <0 
C2.50) 
= K gl ks 86> 0O 


Rewriting the governing equations with the given boundary conditions, 


the equations are 


ag/at = k"(9°6/ar-) (2.34): 
@ = S,0 <0 | 
(2528) 
>o = uf” as dk Gi > 0 
4 aoe §@ = 8) 
(2.36) 
xX > %, (be 5, eco 
25203 Solution 


Considering now equation (2.34) 


ag/ar = k' (a o/ac°) 
writing it in the finite difference form, it becomes 


dg/at] = (gm 6") /ar 
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where pute is one time step ahead of eae 


n 
Initially all the temperatures at the nodal points were assumed 
to be at a constant temperature ue New enthalpies and hence the new 
temperatures, one time step forward, were then found by applying the 
boundary conditions to the governing equations. The time interval At 
was chosen to be so small that only the local temperature and the tempera- 


ture of the adjacent reference points need to be considered in calculating 


the local temperature. It was found that if At was greater than 
i 
2 


explicit computation scheme. 


(ame numerical instabilities resulted because of the use of an 


Briefly the steps are: 


EA Set all temperatures equal to 6 (ie: cbs = 0) and then change 


it into ¢ using equation (2.28). 


¢ =S Oat: stes) =? 0) 
mSia 0 if @ <0 

ii. Apply the boundary conditions equation (2.36) and solve for 
the enthalpy changes at each nodal point using equation (2.39). 
Then using enthalpies from the last time step and these enthalpy 
changes, new values of the enthalpies at each nodal are 
calculated. Enthalpy changes only apply for node 2 to node 
N since this is an imposed surface temperature problem. The 
temperature at the surface therefore is always equal to the 


imposed temperature, which in this case is the air temperature. 
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So for node 2 to (N-1), using equation (2.39) 
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For node N, assuming the heat flux from deep in the ground to 


be negligible 
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iv. Proceed as before, as the time step advances, but now the 

initial temperatures and enthalpies are those just computed, 
sp + 

eu and oo s respectively. 

n n 

As for the surface temperature, 85 = A: 

We The computation is over when the average temperature distribu- 
tion over a year approaches steady state. This numerical 
approach provides onewith a continuous computation of the 
temperature distribution and also the travel of the fusion 
front with a minimum of complications. 

Didnt For the case when the phase change occursover a temperature 
range instead of at a single temperature, equation 234) 
becomes 
; ee i2 
dd/ot = k'(9 6/9) (2.40) 

ime ‘ S 
where k (l-unf/M,.) + (kK g/kKe) unf/M., 6n<20 
(2241) 
a O20 
Kig/Ks = 
and equation (2.28) becomes 
= < 0 
) 5-8 + unf/M. 8 
(2.42) 
Sn 6" £0nh ad Oe "0 
uf S 
where unf is the unfrozen fraction and ue is the soil moisture content. 


aa 


m+1 : 
0 are then the new temperatures one time step ahead of 


n 
ee Both the new temperatures ‘ae and enthalpies fae 


are stored. 


Proceed as before, using new equations (2.40), (2.41) and (2.42). 
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26. 
Jaga fi8) Comments 

In this thesis the steady periodic solution will only be of 
interest, the transient behavior produced by starting from a constant 
ground temperature will not be considered. Rigorousiy speaking, the 
steady periodic solution is obtained ae es the temperatures at all 
times during the yearly cycle are the same in two consecutive years. 
However, since this work is primarily concerned with the temperatures 
averaged over a yearly cycle, a somewhat less rigorous test was used to 
determine if the steady periodic state had been reached. The criterion 
used was that the change of the averaged temperatures from one year to 
the next must be small. Im practice the calculations were always 
ented out through five yearly cycles of temperature variation. If the 
change in the averaged temperatures was greater than 0.02°C between the 
fourth and the fifth years, the calculations were extended to the tenth 
or the fifteenth years as required. 

Input data for this model includes the ground surface 
temperature variation during the year and the physical properties of 
the soil. The ground surface temperature is assumed to be equal to 
the air temperature, which is approximated by the harmonic function 
fitted to meteoralogical data, Appendix I. Temperature conditions of 
Edmonton, Norman Weils and Inuvik were used as the basic cases. For 
soil properties, data compiled by Lunardini (28) was used, Appendix II. 
The basic conditions used for calculating effects of soil property 
variations were taken as a fine grained soil, dry density 1.2 x 10° kg./ 


ae with moisture contents of 15%, 30% and 452. 
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2.3 Results of the numerical calculations 
Dil ed. Temperature distributions in the ground 


Figures (4) and (5) show the temperature distribution in the 
ground for two cases, Edmonton and Norman Wells, both with 30% moisture 
content in a fine grained soil. The temperature profiles were drawn 
for the various time intervalsof the fourth year. 

The temperature profiles for the two cases display the characteris-— 
tic features of all ground temperature profiles, that is the amplitude 
of the temperature variation decreases with depth and the phase of the 
temperature variation lags behind the phase at the surface. The overall 
temperature profiles do not appear to resemble those for the idealized 
analytic model, where Stefan number was zero. It will be noted however, 
for both the cases, the temperature profiles that cross the 0°C line 
tend to be linear from that point to the surface temperature value. This 
is particularly striking in the case for Edmonton where the average air 
temperature of 2.56°C, is very close to the freezing point. Since it is 
the region where phase change occurs, it is thus the most important region 
in Acteenintne the thermal behavior of the ground. Thus, it is reason- 
able to expect that the analytical model will give satisfactory results 
for some calculations. Also, note that in the temperature profiles in 
Figures (4) and (5), the average temperature deep in the ground appears 
to be different from the average temperature at the surface. This is 
the result predicted by the analytical model. 

Figure (6) shows the variation of the average soil temperature 


with depth. It will be observed that the variation occurs only in the 
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ai. 
top few meters of the soil. This depth will be observed later to 
correspond to the depth of penetration of the phase change. This 
indicates that, the phase change or rather the conductivity change 
that results from phase change is responsible for the change in the 
average ground temperature with depth. 

The absolute determination of the ground temperature in 
practice is almost impossible because of the large number of variabilities 
involved. However it is possible to predict the change in the aed 
temperatures that would result from a change in one of its controlling 
parameters. This section will deal with the effects the main parameters 
namely, average surface temperature re volumetric heat capacities of 
the soil Le» thermal heat conductivities ratio kiglke and the soil 
moisture content Moo? have on the differences between the average 


surface temperature cx and the average ground temperature fee 


o 


The'effects that changes in the average surface temperature 


have on the temperature difference a = TS» is shown in Figure 7. Here 
the temperature differences have been normalized by AT, the amplitude 

of the temperature oscillation at the surface. These results are for 

a ratio of conductivities eteen the unfrozen and frozen soil, kuglke> 

of 0.647 which corresponds to a fine grained soil with 30% moisture content. 
It will be noted that for this ratio of conductivities, the temperature 
below the active layer is always colder than the ee tees surface 
temperature and the maximum difference occurs when Bs is slightly 


greater than the freezing temperature. The theoretical curves indicated 


by the solid lines in Figure 7 meet at a value of (T-T {AT = —0.14. 
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The curve on the left corresponds to a soil temperature deep in the 
ground that is above 0°C and the one on the right corresponds to a 

soil temperature below 0°C (that is permafrost grounds). The intersection 
of these curves and thus the minimum value of the parameter (T,-T,)/4T, 
corresponds to the condition under which the greatest difference exists 
between the average surface and ground temperatures. The minimum value 
of (To-T,)/at for this case is -0.14, which for a typical temperature 
range AT = 20°C would correspond to a difference between average surface 
and ground temperatures of 2.8°C. The position of this minimum value 
and thus the minimum value of (Te-T,)/aT is a complex function of the 
conductivity ratio ke Ke: It can be calculated by equating equations 
(2.22) and (2.24) of the analytic model. 

The results of the numerical calculations for Le a WA Aa 10°, 
that is mM. = 30%, agree with the analytic results except near the 
minimum point in Figure 7. Here calculations must be run through 15 
yearly cycles before a steady periodic state is reached. 

An increase in the unfrozen moisture causes a decrease in 
(1 -T,)/ OT, as indicated in Figure 7. This is due to a decrease in the 
range of conductivity. 

A change in the moisture content of the soil effects both the 
latent heat of fusion Ly and the ratio of the conductivities K Ke 
in the soil. These effects of the two parameters will first be looked 
at separately. 

Figure (8) indicates that if the calculations are carried out 
through enough yearly cycles, approximately 15, the value of (T,-T,)/aT 


obtained is independent of the volumetric heat capacity of the soil Le. 
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The result is in agreement with the analytic model. 


The ratio of conductivities, k ./k,, has a much more 


Chee wea cy 


significant effect. The relationship between Ke! Ke and (T,-T,)AT 
is linear, that is as k gl Ke increases, (T,-T,)/aT increases, as 
indicated by Figure (9). For every unit ee ae in k gl ke» there 
result a corresponding 0.166 units increase in (1 -T,)/OT. This. 2s 
at a value of (Tete) / aT = 0.341. 

The graph of (T,-T,)/AT versus the moisture content Moo in 
Figure (10) will only confirm the results of the two previous findings. 
The graph shows that an increase in moisture content will result in an 
increase in the magnitude of (T,-T,)/AT. This results because of the 
increased difference between kit and ke that occurs at the higher 
moisuture content. Also changes in moisture content have less effect 
on the ground temperature in the colder region such as Inuvik than the 
Sree ale Goines Pee onwench as Edmonton. This is because as shown 
in Figure (7), the Brrects bf conductivity differences between frozen 


and unfrozen soil are the greatest when the soil is near 0°C. 
2.4 Depth of the active layer 
DOL Depth of frost or thaw penetration 


Figure (11) shows the frost or thaw penetration depth during 
the fourth year for the three locations, Edmonton, Norman Wells and 


Inuvik, with 30% moisture content. The method of obtaining the frost 
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or thaw penetration depth will be discussed briefly below. 
Dy heal: Method of calculating the depth of thaw or frost penetration 


Consider first the case where there is frost penetration. 


That is a frozen phase over an unfrozen phase. 
ig Te 703 then i 18) 
n n-1 
and the depth to the fusion front is 


a = 
y n OFD oF 


The 0.5 occurs because the first nodal point represent only 1/2(Ax). 
The latent heat of the nth nodal point is - and it represents the 
fraction of that element unfrozen. For thaw penetration, that is an 


unfrozen phase over a frozen phase. 
Eat el =e } then T a0) 
and the depth to the fusion front is 
ys =n- 1.5 + om 
The depth of thaw or frost penetration, y, is given by 


= y* 
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Example: 
Considering now the Edmonton case, that is frozen/unfrozen. 


Time = 4th year 


Nodal Points Temperatures - Enthalpies 
Ground 


i. Surface 
-12.89848 =U oS 


°2 =10'. 79542 -0.14394 
Je =O. f 4095 =O7LLGS5 
“4 ae 1 Gia) 4 =0'. 09033 
= ie 
“5 Cay = Ano Lig -0.06562 
—s Fur AG 
°6 = LOO0I8 -0.04255 
f, =. 56493 -0.02087 
ns) 0.0 0.94294 
_-—— eS eared en ¢ = 1.000 —- — —— — 
°9 Unfrozen 0.18144 1.00363 
Region 
¥eu="Se-00.5 — 0.94294-— 
= 6.56 
Veer Fe Oe 
At = 0.06 
am we ke E AT 
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Hence the depth of frost penetration, 


(6256s 0-06 x 2.6713) meters 


<< 
It 


1.05 meters 


As indicated in Figure (8), for the stated conditions, that 
is imposed surface temperature, average temperature TN = 2/556" C, 
amplitude of oscillation AT = 15.96°C and with 30% moisture content, 
the maximum frost penetration depth is 1.8 meters, about 5.9 feet. The 
generally accepted and known frost penetration depth in Edmonton is 
about 8 feet. For the cases of Norman Wells and Inuvik, both of which 
are naturally permafrost locations, the maximum depths of thaw are 1.225 m 
(about 4 feet) and 1.00 m (about 3 feet) respectively. It will be 
observed that these depths of the active layer correspond to the depths 
over which the variation in the average ground temperature occurred in 
Figure (6). This is the expected result, infthat it is the conductivity 
variations that accompany the phase change that causes the average 


ground temperature to vary with depth. 
Detee Effects the various parameters have on the active layer 


As it is possible to predict the change in the ground 
temperatures, it is also possible to predict the change in the depth 
of the active layer that would result from a pHaaee in one of its 
controlling parameters. 

Figure (12) shows the effects the changes in the average 
surface temperature have on the depth of the active layer. A collation 


with Figure (7) shows that they exhibit the same behavior with the 
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maximum depth occurring when the value of Cea] AE ==O.V14. “For a 
conductivity ratio k gl ks of 0.647, which corresponds to a fine 
grained soil with 30% moisture content, the maximum depth of the active 
layer is 1.85 meters about 6.1 feet. 

Again the effects of the latent Rea of fusion and the 
conductivities of the frozen and the unfrozen soil will be examined 
separately and then in combination as they are related to the moisture 
content of the soil. 

For both the naturally permafrost location and the non 
permafrost location, the depth of thaw and frost penetration respectively 
decreases as the latent heat of fusion Le increases, as indicated in 
Figure (13). For the Edmonton conditions, case B, the decrease in active 
layer depth closely follows the theoretical predictions of the analytic 
model. The decrease in depth is therefore proportional to (L_) ie 
For the Oe Fite permafrost conditions, case A, the numerically predicted 
active layer depth fottows tHe same trend as predicted by the analytic 
model, however the actual values from the numerical calculations are 
20 to 25% less than those predicted analytically. The accuracy of the 
analytical prediction is in general best, as shown in Figure (12), when 
the average ground temperature is near the freezing point. This is 
because in this condition the Stefan number is near zero - the condition 
for which the analytical model was developed. 

The analytic model would predict that the depth of the active 
layer for the permafrost soil is dependent only on the conductivity of 
the unfrozen soil and that for the non-permafrost soils the depth would 


depend only on the conductivity of the frozen soil. In particular the 
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LATENT HEAT OF FUSION, Ly {x108 joules/m? } 
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FIGURE 13 — The effects the latent heat of fusion L¢ has on the depth 
of the active layer, for iigKe = 0.6471 


a 


“a 


” 


eth 


ee pele ny 
jp naan) i 


sig _ ot 
Rat ae 


4 
yp ine ie 
Mee 
=) 
> 
Nea, © 
. 
yo 
& 
Teo “ 


~ 
34 
* 
. 
‘ 
i 5 
— 
*, ‘ 
v4 
At es 


- 
+o 

) 

4 Vi 
ae J 
« e 
~ 


a4 ae 
iA © 
‘Fe “ ce 


+} aie Bie wv 


ak 

ty f 

ay he ees 
ie ae ers 


OP < 


Tete 


Bare 


Linde 
a 
: ‘ 
age 
afi Oh, abe rhe 
j 
j 
‘ a) mpl 5 
s 
i 
: 
te 
s¥ 
e +o Les 
2 
~ rs F 
Awl = 
.< 


He vial 9 a nn dint ay “ “ing event 
pe) oy aml " (ah 


Weg et steel ‘ Ne 
Ve a f ; 
4 - “he bi , he ey 
- , , t 
~ fe ‘ if y * ’ MA 
| 7 iy ie 
AS ; ‘Ge 2 F a, 
y * =} 


45. 


depth in each case should be proportional to the square root of the 
appropriate conductivity. Figure (14) shows the results of the 
analytic and numerical models. 

Figure (5) shows the relationship between the moisture content 
of the soil Moo and the depth of the active layer for the various 
locations. For the naturally permafrost locations such as Inuvik and 
Norman Wells, the depth of the active layer decreases as the moisture 
content increases. As for the Edmonton case (non-naturally permafrost 
it is just the opposite, that is the depth of the active layer increases 
as the moisture content increases. 

The results for Edmonton are interesting, in that just from 
a consideration of the latent heat which increases with moisture content 
one would expect the depth of the active layer to decrease. The increase 
in moisture content however also increases the conductivity of the 
frozen ground and this acts to produce the increase in the active layer 
depth observed. For the two permafrost locations however, the effect 


of the latent heat dominates the behavior of the active layer depth. 


That is as the latent heat increases the active layer depth decreases. 


2.5 Conclusions 


The analytic model developed, based on S- = 0, gives results 
that are in close agreement with the numerically calculated results. 
The results obtained demonstrate that there exists a difference between 
the average ground surface temperature and the average temperature just 
below the active layer and that this difference is due to conductivity 


changes associated with phase changes. The difference between the two 
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FIGURE 15 — Variation of M,, with depth of active layer 
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mentioned temperatures is largest for soil with large moisture Sentene 
and in geographical locations where the average surface temperature is 
near freezing. 

Interestingly enough, both the analytical and the numerical 
models show that the average ground temperature below the active layer 
is independent of the latent heat of fusion. The depth of the active 
layer is inversely proportional to the een e heat to the one-half power 


as would be predicted from the classical Stefan problem. 
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PART III 


WITH AN INSULATING LAYER OF VARYING U-FACTOR 
Se LierToductLon 


Notnally, instead of the imposed surface temperature situation 
discussed in the last section the actual ground surface has an insulating 
cover (vegetation cover in the summer and snow in the winter). This 
insulating cover will have a major effect on the ground temperature. 

This section of this thesis will examine the effects of this insulating 
cover on the ground temperature. 

The analytical model presented will neglect the effects of the 
latent heat of fusion on the average ground temperature and will ewes 
on the effects of the insulating layer. In the numerical model all the 
factors will be included and then a collation will be made between the 


two models. 


3.2 Formulation of a ground temperature model with an insulating layer 


of varying U-factor 
Seal Analytical Model 


Assume that the infinite ground mass can be replaced by a plate 
of some effective heat capacity pc. This plate is then assumed to be 
at the temperature of the ground surface and that no heat flows into 
or out of it from the bottom. The effective heat capacity of the plate 
would be chosen so that the plate will have the same thermal response 
as the ground normally would. Choosing an actual value for pc, will 


be discussed later. This approximation however neglects the storage 
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of energy in latent heat which will affect the results as will be 
observed. 
The model consists of a layer of insulation overlaying, a 


plate of finite heat—capacity, Figure (16). 


at 
FIGURE 16 
T, = surface temperature 
= w TPodAh Gaiyipe 
U(t) = U-factor of the insulating cover 


k/d 


U(t) will be assumed to be of the Coen U + AU sin(wt + o>: 

The phase angle ge is the difference in phase between the varia-~ 
tion of the U-factor and the surface temperature. A value of oF = 0 means 
that U(t) is in phase with the surface temperature, that is the snow cover 
in winter is a better insulator than the vegetation cover in the summer. 


The normal case one would expect would be o near zero. When oe = T 
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U(t) is 180° out of phase with the surface temperature, implying 

that the summer vegetation cover is a better insulator than the snow 

in winter. This case would not be as likely, however it could possibly 
occur if the ground were covered with a thick layer of moss which became 


very dry during the summer months. 
ae ae Solutions 


For this model the governing equation for the ground surface 


temperature eae is 
pe, (dT /dt) = U(t)[T, (t) = cae (371) 
where T(t) = vs + AT sinwt G22) 
U(t) = U + AU sin(wt+9 ) Ges} 


and T isgsaeunction «of time, :t. 


gl 
Let 
0 = Cra) = T(t) 
° d0/dt = dT dt - dT /dt 
hence dT /dt = de/dt + dT /dt 
but dT. /dt é wAT coswt 


which means that 


dT _/dt = d6o/dt + wAT_ cosut (3.4) 
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Substituting equation (3.4) into equation (3. 


pc. (de/dt + wT coswt) = - U(t) 86 


or 


dé/dt + U(t)9/pc, = - wAT coswt 


from equation (3.3) 


U(t) U + AU sin (wt + oD 


= U + AU sinot cos | + AU cosuwt sing | 


let 
iT = 
U/oc, A 


tt 
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AUcos¢ /oc, 


I 
(@) 


AUsing /pc, 
Equation (3.5) will now become 
d6/dt + (A + Bsinuwt+ Ccoswt) 9 = - wAT coswt 


Assume 


6=a +. {fa sin(awt) + b_cos(mut)} 
ayn n 
n= 


Sa b/dt = =. {a nwcos (nwt) - b nwsin (not) 
n=1 
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Substituting equations (3.9) and (3.10) into equation (3.8), the- 
expression becomes 


co ao 
X {a _nweos(nwt) - b_nwsin(nwt)} + [a + 2 {a sin(nut) 
Tels. tl n Oe oat 


+ b cos (nut) }] x {A + Bsinwt + Ccoswt} = -wATcoswt (3.11) 


Multiplying out the second term on the right hand side and then regrouping 


the common terms, equation (3.11) becomes 


Cc 
Aa + Ba sinwt + Ca coswt + 2 {(Aa -nwb ) sin(n)wt 
Oo re) Oo af n n 


+ (Ab_tnwa _) cos(n)wt + 1/2(a_C-b_B) sin(n-1)wt 
+ 1/2(a_Btb, C) cos (n-1l)wt + 1/2(b_Bta C) sin(ntl)wt 


+ 1/2(b_C-a_B) cos (nt+l)wt = -wAT coswt (3.12) 


or simplifying further by letting 


Constant term 
= 0 
Aa, + 1/2 (Ba, +Cb,) 
sinwt term 


cz = 0 
Ba, + (Aa,-ub,) + 1/2 (Ca, Bb, ) 
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cosw term 

Ca, + (Ab, twa, ) + 1/2 (Ba,+Cb,) = -wAT 
sin 2wt term 

1/2 (Bb, +Ca,) + (Aa,-2ub,) + 1/2 (Ca,-Bb,) = 0 


cos 2ut term 


- me + - < 
1/2 (Cb, Ba,) + (Ab, +2wa,) | 1/2 (Ba, cb.) 0 
sin(N)wt term 
1/2(Bb,,_,+Ca,_,) + (Aay-Nub,) + 12 Cae Bb a) 0 
cos(N)wt term 
-Ba vy ( y = m. 
1/2(Cby_ 5 B ae + (Ab. Nua.) + 1/2 (Ba, +Cby 14) 0 (3213) 


Writing these in the matrix form 


[ A B/2 C/2------------ 0 
B A -w C/2 -B/2 ------ 0 
C w A B/2 C/2 ------ 0 
0 C/2 B/2 A  -2w C/2 -B/2 -0 
0 -B/2 C/2 2u A B/2 C/2 -0 
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The matrix in equation (3.14) is of infinite dimensions and 
thus exact values of the coefficients cannot be obtained in this way. 
However, under certain conditions the matrix may be approximated by 
truncating it at a certain value of N. In particular it was found by 
expansion of the first few terms, that Ker U/ec.w and AU/U are 


small parameters, the coefficients of 'a' and 'b' become small for large 


values of N. The approximation obtained by setting ay and ens 0) 


for n > 2 will be obtained here. The matrix will be reduced to 


A B/2 C/2 | ay 0 
B A -W | ay = 0 (3.15) 
C Ww A 1 -wAT 
Solving for ay» a, and by 
fio B/2 c/2 
0 A -W 
-wAT W A | 
Bone A B/2 C/2 
B A- =W 
Cc Ww A 


Mi -wAT(-wB/2-CA/2) 
0 A(A*4+u2) - B/2(ABHCw) + C/2(Bu-CA) 


AT B/2A (1 4 CA/Bo — 1/u*(AC-B-/2 & c7/2)} 
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hence 
hat = sing —) 9) 
eee cost. ie ee 
0 u pc _w cos 2 2 
(ec w) 2(pc_w) 
e e 
(3.26) 
Similarly 
A 0 Gi 
B 0 -W 
a G -wAT A ie 
“al A B/2 C/2 (3.17) 
B A —W 
Gc W A 
2 Cc 
7 A(- w AT) + 2 (- BwAT) 
aeRO Fan pier ete Corwin Ton ks 
A(A'+w) - > (AB+Cw) + 3 (Bw-CA) 
2 2 
‘ ACEI, $2) da an 
Soult las ai 5 (A 5 ae 
Ww 
hence 
a a = Wy 
I 
a, =~ AT + AT aoe =—— cus 
p a pc w eCn U 
Det 
por) Sil 20 ce} (3.18) 
en ca u 
U 
and finally 
A B/2 0 
B A 0 
Cc W -wAT | 
= Shits ah ent 
eT A B/2 c/2 (3.19) 
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Neglecting all the higher terms, ie. terms with U/oc,w)” or higher 
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U Tee er Ane 
b = AT Gel {(-1) + 4 = + 4) cos 26 3 (3.20) 


From equation (3.9) 


foe] 
6 = a+ f {a_sin (nwt) + b cos (nut) } 


+ ths = 
0 + a, sinwt 4 by COS UME Vet o's rte «cus cee ee eees (3221) 
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Substituting equations (3.16), (3.18) and ©.20) into equation (3.21), 


the following expression follows: 
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To obtain the difference between the average temperature (ToT) 


equation (3.22) may be integrated over a time cycle. The sinwt term 


and the coswt term will drop out, leaving the following expression: 


but 


Qn 
f° @ dt = 6(21/w) 
oO 
hence 
i cos¢ = sing ee 72 
@ = AU AT 5 Ya —— Te - + } 3G723) 
U Cee u (pc wv) 2(pc.w) 


In the limit U/ec.w = 0, noting that AU must by definition be less 


than U, equation (3.23) gives 
T i-T = U 2 3.24 
a To )/AT (AU/U) (cos¢ / ) ( ) 


Also note that from equation (3.22), that this limiting situation implies 
that the ground temperature variations with time are negligible. in 
this case the temperature difference Tat, varies linearly with 


AU/U and as the cosine of the phase angle oe This implies that 
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for oa = 0, that is for a snow cover of good insulating value, the 


average ground surface temperature ca is greater than the average 


surface temperature To: For o = 180°, that is the summer vegetation 
cover is a better insulator, the average ground surface temperature 


cea is less than the average surface temperature Tt There is no 


difference between T and T if however @ = 90° If U/oc Qs 
gl Ss u e 
not zero, that is if the fluctuations of the ground temperature are 
accounted for, equation (3.23) shows that a difference between To 
and i exists at 8 = 90°. The amount of the difference is 
Tiare? = (AUxAT)/2pc, w (8525) 


For larger values of U/ec.w and AU/U, higher order terms become more 


‘significant and equation (3.24) is not a good approximation. 


3.3 Numerical Model 
The method adopted here will be similar to that developed for 
the imposed surface temperature except for a slight modification at the 


soil surface level. 
S Pee re The difference between this and the previous model 
(a) At the ground surface with imposed surface temperature 
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FIGURE 17 (a) 
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FIGURE 17 (b) 


At the soil surface 
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(3.26) 
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one 
Ree kpt AT/peL 


« FUS = Ux /k, (3.28) 


Initially the ground temperature were all assumed to be equalled to sis 


or aT = 8, ae teleners os = 9. 
Enthalpy 
Ad, = {FUS (6 .-6 , )2AcxRIF - (8-8) 2xRT} (83239) 
Z 
where RTF = At/At (3.30) 
2 ; 
and RE. =) At/Ac if 8) he 8, <0 
(3.34) 
* 2 , 
= (kK g/Ke) At/At tf FT + 8, a0 
ire em m 
oy : 9) a ma 
Temperature 
ele tL iia 
6, = oy /S_ .f, 4 <a) 
(3432) 
eo nt. m+1 
GR Se if $) > 0 


The main difference between the two problems is, instead of 
setting the ground surface temperature equal to the air temperature, the 
ground surface temperature will now be varying. This is caused by the 


introduction of a varying insulating layer which separates the atmosphere 
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from the ground surface. The temperature on the surface of the insulating 
layer is now assumed to be equalled to the air temperature. 

So, except for these modifications at the ground surface 
level, all the proceduces and parameters remain unchanged. The varying 
insulating layer is approximated by a harmonic function which may or 
may not be in phase with the air temperatures. Again temperature 
conditions of Edmonton, Inuvik and Norman Wells were used as the basic 


cases. 
3.4 Results of the Numerical Calculations 
eal! Effects the insulating layer haveon the various temperatures. 


The temperature profile, as illustrated in Figure 18, display 
similar characteristics as those for the imposed surface temperature, 
except that the insulating layer decreases the amplitude of the 
temperature variation at the ground surface. The amplitude of surface 
and ground surface temperature variations for the three cases, together 
with the soil and insulating properties are listed in Table 1. The 
properties of the insulating layer chosen for the three cases are rather 
arbitrary as there is very (ieee data available on the insulating 
values of the various surface covers. Some rational range of values 
used is given in Appendix III. In all cases the phase angle oe Lee0. 0; 


that is the situation in which the winter snow cover is better insulator 


than the summer vegetation cover. 
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TABLE 1. Amplitude of Temperature Variations 


Locations, soil and Amplitude of Amplitude of ground 


Insulating properties surface temperature surface temperature 


variations (°C) Vabiatzons =¢ C) 


Edmonton: 


T. = 2.56°C, AT = 15.96°C 
Ue 250, AU/U = 0.6 18.1°C: fa =1218°C 154°C REG 503°C 
¢ = 0.0, M_ = 30% 

u so 


Normal Wells: 


Ts = =O 22°C se 22 106°C 

Wi (0.75, AU/U. = 0.8 15 4li2@ to: =27%6°C Glee £0) 26596 
$, 7 0-0, Mo, = 30% 

Inuvik: 

us = OCOh Cu, = 22.042 °C. 


0.75, AU/U = 0.8 12. 0° Catos—-3ie ce AeO Ge foe 10 27 


The amplitude of temperature variation at Norman Wells decreases 
the most because its average ground temperature is closer to OC, the 
temperature where the phase changes occur. As indicated by Figure Choy, 
the temperature deep in the ground is greater than the average surface 
temperature T.: This is true for all the three cases. Hence with 
the introduction of an insulating layer in phase with the air temperature, 
the ground temperatures become warmer. This is the result predicted by 


the analytical model. 
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Figure (19) shows the variation of the average ground 
temperature with depth, with an insulating layer in phase with the 
air temperature. The temperature profiles show that the average 
temperatures deep in the ground are less than the average ground 
surface temperatures. These results are in agreement with and have 
been dealt with in the previous section where the surface temperatures 
were imposed. The temperature profile also showsthat the average 


ground surface temperatures eat are greater than the average surface 


temperatures Tee These results are predicted by the analytical model. 
The differences between the two average temperatures ce and the are 
very Significant for both the Norman Wells and Inuvik cases, both of 
which are in the naturally permafrost locations. Also indicated are 
that the differences between Tig and Ty are very much greater than 


that between Th and Ty implying that the effects of the insulating 
layer are much greater than the effects of the thermal properties of 


the soil. 


Sey Effects the insulating layer have on (To-T, 1) 


This section will deal mainly with the effects the insulating 
layer has on the difference between the average top surface temperature 


T, and the average ground surface temperature Ta The main parameters 
of the insulating layer U are, the average insulating cover U, the 
amplitude of variation AU and the phase angle oe It should be 

noted that the theory used in the figures would only be the "simple 


theory" or equation (3.24), which is the dominant term of the more 


complicated theory or equation (3.23). 
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Equation (3.24) or the "simple theory" (where U/oc.w = 0.0) 
is 

0 = (AU/U) x (cos® /2) 
Also, oe in the analytical model is o x 27 in the numerical model. 
The main reason why the more complicated theory or equation (3.23) was 
not used as the standard theory was because of the difficulty in 
predicting the value of the effective heat capacity of the soil pc. 
This is because the effective heat capacity pc, varies as any one of 
the controlling parameters varies. Attempts were made to obtain a 
suitable average value of pc, and substituting it into the more 
complicated theory to obtain a theoretical curve that would match 
the numerical results more closely. This was not possible in general 
because of the sensitivity of pc, to the other parameters of the problem. 

The.effects the ratio AU/U has on (2,71) /AT are shown in 
Figure (20) for two cases, one when the insulating layer is in phase 
Cee = 0.0), while the other when the insulating layer is 180° out of 
phase (os = 0.5) with the air temperature. For the case when (ae =e 
(Toyo) /OT increases almost linearly as AU/U increases and the 
numerical result is about 45% less than the results of the simple theory, 
that is neglecting the terms in U/pc.u. When oe =. Ori5:J (Toyo Td /aT 
decreases as AU/U increases. The numerical results here are about 
25% more than the results of the simple theory. The big difference 
between the simple theory results and the numerical results is mainly 
because of the assumption in the simple theory that U/ec.u = 0.0. An 


inspection of the more complicated theory, equation (3.23), indicates 
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0.5 | 


THEORETICAL (simple) y, 


0.4 (4,7 0.0) 


0.3 
NUMERICAL 
O (p= 0.0 ) 


0.2 


0.1 


20:1 


NUMERICAL 
-0.3 EDMONTON @ ($,70.5) 
PE2.562°C 
AT=15.96°C THEORETICAL (simple) 
-0.4 30% Meo | (p°0.5) 


U=0.75 watts /m*-°K 


705 


FIGURE 20 — Variation of AU/U with (Tg — T,)/AT 
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that any finite value in U/pcw would decrease the absolute value of 
(ToT) / AT hence closer to the numerical results. 

Figure (21) shows the relationship between (T,477,)/AT and U. 
for a given value of AU/U. This insulating layer model will approach 
that of the imposed Mitace temperature tty when U approaches infinity. 
Analytically, the simple theory indicates that (T4774) /A7 is independent 
of U but however the more complicated theory qualitatively shows a decrease 
in (Tort /at as U increases. The numerical results confirm this 
behavior. They also show that the decrease in (Toy7T,/a7 is greater 
at 15% moisture content than at 30% moisture content. This could be 
attributed to the fact that the effective heat capacity (pc.,) is greater 
in the 30% case. 

Figure (22) shows the relationship between (Too) /aT and 
the insulating layer phase constant tie for the various moisture content 
at Inuvik and Edmonton. When phase constant ¢ = 0.0, the value of 
(Toyo /at is positive se it increases with moisture content. For 
phase constant ¢ = 0.5, the value of (TT) /at is negative and it 
decreases with moisture content.- However, when the phase constant 
ae = 0.25, the simple theory indicates a zero value for (To7T,)/4T 
but the numerical results shows that the value of (T5477) /aT is 
positive and it decreases as the moisture content increases. On the 
other hand the more complicated theoretical solution predicts a positive 
value of (Toyt,/At when o = 0.25. The above observations are true 
for both Inuvik and Edmonton cases. Generally the values of (T,47T,)/ 47 


for Inuvik case are greater than those for Edmonton at all phase constants. 
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~~ 
EDMONTON ™ 
T= 2.00. © ™~ 
ee 
AT =15.96°C ee 
$= 0.0 
AU/U= 0.8 


—————— NUMERICAL RESULTS FOR 30% mep 
——-—— NUMERICAL RESULTS FOR 15% m, 


1.0 20 3.0 40 5.0 


AVERAGE U-FACTOR, U, watts /m?~°K 


FIGURE 21 — Variation of U with (Tz —T,)/AT 
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and that the moisture content effects the values of (T 


ae at 


Inuvik significantly more than those at Edmonton. 
3.5 Depth of the active layer 
S561 Depth of frost or thaw penetration 


Shown in Figure (23) is the typical thaw penetration for 
Norman Wells. Depth of maximum thaw or frost penetration for the three 
locations under consideration are listed in Table 2. Listed also are 
the soil and insulating properties. The method of obtaining the depth 
of frost or thaw penetration is similar to the one discussed earlier 


in Section (2.4.1), for the imposed temperature case. 


TABLE 2. Depth of frost or thaw penetration with 


imposed surface temperature 


Locations, soil and insulating properties Depth of thaw or frost 


penetration (meters) 


Edmonton: T. = 2956-°C, ‘AT = 205.9646 
U = 2.0, AU/U = 0.6 0.925 m. (= 3 feet) 
ie 0.40; Mo 1570 

Norman Wells: c =y-6.22, AT = 2206-6 
UF=.0<75; AU/U ="088 0.825 m. (= 2.7 feet) 
oe= 020;.M.. = 302 
u so 

Inuvik: te = =O061°C, ATi=s 2200250 
U = 0.75, AU/U = 0.8 0.425 m. (= 1.4 feet) 
¢ = 0.0, M. = 30% 
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Snide Effects of the insulating layer on the active layer: 


Figures (24) and (25) show the relationship between the 
average insulating layer U and the depth of active layer for Edmonton 
case, with various phase constant On and: moisture content Moo" 
Generally, if all the parameters remain unchanged, an increase in U 
will cause an increase in the depth of the active layer for the Edmonton 
case of the non-naturally permafrost location as indicated by both the 
figures. Specifically for a given U and moisture content Moo? the 
depth of the active layer increases as the phase constant oF increases 
and for a given U_ and phase constant op the depth increases as the 
moisture content decreases. 

For Edmonton with 30% moisture content, U = 0.75 and phase 
constant op = 0.0, the depth of the active layer decreases as AU/U 
increases as indicated by Figure (26). However, if the phase constant 
ta = 0.5, the depth of the active layer increases at first gradually 
and then rapidly as AU/U_ increases. 

Figure (27) shows the relationship between the moisture content 
Moo and the depth of the Spee: for both Inuvik and Edmonton at 
the various phase content dee Figure (27) indicates that for the Edmonton 
case (non-naturally permafrost location), at a given U-factor, the depth 
of active layer increases as as increases and decreases as the 
moisture content ee increases, the result obtained in Figures (24) 
and (25). For Inuvik, a naturally permafrost location, for a given 
U-factor and phase constant Oe the depth of the active layer decreases 
as the moisture content increases, similar to that for Edmonton but 


the depth decreases as the phase constant hy increases, opposite to 
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FIGURE 24 — Variation of U and ¢y with depth of active layer 
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AVERAGE U-FACTOR, U, watts /m=°K 
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FIGURE 25 — Variation of U and M.9 with depth of active layer 
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FIGURE 26 — Variation of AU/U with depth of active layer 
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SOIL MOISTURE CONTENT, m,,,°% 
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FIGURE 27 — Variation of M,, and ¢,, with depth of active layer 
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that for Edmonton. 
3.6 Conclusions 


With the addition of an insulating layer, the amplitude of 
temperature variation at the ground surface is smaller, the greatest 
change occurs when the average ground temperature is closest to 0°C. 
Except for the smaller amplitude of temperature variation, the profiles 
of the ground temperature are the same as those for the imposed ground 
surface temperature case. An insulating layer in phase with the air 
temperature would cause the average ground temperature to be warmer 
than the average air temperature. If U increases, that is approaching 
the imposed ground surface temperature condition, the effect of the 
insulating layer on the ground temperature decreases. The effect of 
the insulating layer is largest for high moisture content soils. 

For both Edmonton (non-permafrost location) and Inuvik 
(naturally permafrost location) cases, an increase in moisture 
content Moo would cause the active layer depth to decrease. However, 
for a given Moo? an increase an oe would cause the active layer 


depth in Edmonton to increase whereas that in Inuvik to decrease. 


4 i Ne ; 
e j 1 
~~ 
' . ES Se 
me oh ae ie 
rai Sk.@ COE EP tb 
Hedy oko. anaeegie See 
ca > 1 
i _ ; snk its iat sina | ey hae rH 
wid bh ira mt tal m: FOWL nah soleil 
‘diay wi Be a idee panty ae 
- eu z _ +) pee NN iis 
(qi ot aka Akeamesik, Ay’ at pale ves 
pe 4 
iz pwGPs mt? inh aes aanohnr, saan is 
tie itt yin py spleen verre rn ie 
2 ' bw ‘ fi y A 
he ‘a Aga | Hebe = ila Fe shagnph dl 
Cy py, To 
yah 0 atv tisuastl sec taaaa sega en, 
y ~ ‘ ay oe 7 a 
pha sani a ace nat ‘ners a ‘Sos hebea (uh i woes 
ee ee i a Ae, : 
yay pier e (edits cine nim ‘wt ‘ovine been cae ‘gi 
+ Ee eh pa : yA 
pais mielged, aA’ Tid Whisd i: a) a, ‘owabrun oy. Mygee pate on 
2G " ca 2S wits ‘ he Seine dasesosd a8 pitas ai. d 
* Was a 4 Ri a Was oy mn a : 
Aw 7h i , ; ro ' 7 a 
F my ; . ee / ; = ' a Pal 
ey ) ) 7 
r ' , j wd eP ; . 7 
} | ) a ~ . a 
t 1 oat Ay a re a : 
a a a 
‘ . fi a ee : : : 7 ; bh un 
y ; i 4 : Ma ¥: 
L a P ' nS ‘ 
: =k Vi : 4 ae ; 
4 ma hg : 
: “! a i‘ . rs | ae ae wy a 
- : nae ne an j 4 =i. = ov 
4 . ' “ - x ni js } - r ve ¥ ; 7 A 
fs = Fy reeds “ip ers Fi 7 i 
= La ; “ a! i a 
: Len a] eS; Aa ; van ae 
- ay Taree et , bs 7 pear i _ - 


PART IV 


UNSTEADY HEAT CONDUCTION WITH CHANGE 
OF PHASE AND VARYING HEAT FLUX AT ONE 


BOUNDARY 


4.1 INTRODUCTION 


In this section, instead of assuming Lee equals the air 
temperature To Ts here will be calculated from a balance of heat 
fluxes at the surface. It is reasonable to assume that the calculated 
T, will be different from the air temperature T. because of the 


various complex heat fluxes. 


Je a Atmosphere 


Insulating 


FIGURE 28 


Figure (28) shows the various major average temperatures in 
the ground-atmosphere system, 
briefly, Lee # To : - because of the freezing and thawing 

Th ? cee : - because of the insulating layer 


- £7 : - because of the heat fluxes 
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4.2 Disposition of Solar Radiation 


Surface 


FIGURE 29 


Figure (29) shows the major components of the disposition of 
Soler radiation in the earth-atmosphere system. The components are: 
i) Ro» incoming solar radiation 
ii) rR» amount of Ry reflected 
ii1.) Ry long-wave radiation down from the sky 
iv) Ry? long-wave radiation (infrared radiation) from the ground 
v) S, sensible heat flux 
vi) a latent heat flux 
vii) Co» heat flux conducted through the insulating surface into 
the ground. 


All of these components will be individually discussed later in this 


section. 
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The air temperature T_ will be assumed to vary as a pure 


a 


harmonic function of time around an average value (Appendix I). Consider 


only the case for Norman Wells, the average air temperature T. is 


-6.22°C and its amplitude of oscillation AT.» ise27.06C. 
i = qT. At sin(wt+9,,) 


The differentiation between summer and winter conditions will be 
determined by whether the air temperature qo is above or below the 
freezing temperature 0°C. With reference to Appendix I, for Norman 


Wells 
when Oe4 2st 083 (summer conditions) 
0.42 2 tn > 0.83 (winter conditions) 


Heat Fluxes 
Sea ke Incoming Solar Radiation (Ro) 


This value of Ro is usually the measured value and is 
independent of the nature of the ground surface. The incoming solar 
radiation is made up of diffuse and direct radiation and is the 
total solar radiation that reaches See earth surface. On the average 
about half of the solar radiation intercepted by the earth eventually 
reaches the surface. The rest is either reflected and scattered back 
to space by clouds and atmospheric constituents or is absorbed by clouds 
and atmospheric constituents. Appendix IV shows the range of short wave 


radiation (R_) for a latitude of 65°N. This data was recorded from 1943- 
Ss * 


1957. 
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These values are assumed to be those of Norman Wells (65° 17' N) and 


that the terms are assumed to vary sinusoidally about an average value. 
=R + + : 
R. AR. sin(wt op) (4.1) 


US Are Reflected Short-Wave Radiation (rR,) 


All solar radiation incident on the earth surfaces is not 
absorbed there, a certain portion of that is reflected and is lost to 
Space. The amount reflected depends primarily on the color, moisture 
content, texture and the composition of the surface. In general, a 
wet or dark surface has a lower reflectivity than the dry or light 


surface. 
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TABLE 3. Albedo of Natural Surfaces (r) 


ners 


Type of Surface Reflectivity, xr 
Fresh snow cover Ook US is) 
Dense cloud cover 0260 = "0.290 
Qld snow cover 0.40 - 0.70 
Clean firm snow 0750 =" 0.65 
Light sand dunes, surf 0.30 = 0.60 
Clean glacier ice 0.30 - 0.46 
Dirty firm snow On 20r= 0.50 
Dirty glacier ice 0.20 - 0.30 
Sandy soil O.15> = 0540 
Meadow and fields OeL2=n0.50 
Densely built-up areas 0.15 = 0.25 
Woods 0705 = 0.20 
Dark cultivated soil 0.07 =. 0220 
Water surface, sea 0.03 — 0.10 


The values given in Table 3 are obtained from Geiger (45). 
Observing the values given in the table, shows that the albedo of 
different natural surfaces vary from about 3 to 30%. This is because 
the albedo varies with the moisture content of the soil and also varies 
both with the wavelength and the angle of incidence of the solar rays. 
The values of the last two cases are given in Sellers (43). In the 
case of cultivated soil, albedo depends strongly on the moisture content 


of the soil, an increase in soil moisture decreases the albedo. For 
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the standard data set used in the calculations, the reflectivity was 
assumed to have an average winter and summer values, r and ts: 


These values were chosen to be representative of vegetation in summer 


and snow cover in winter, tae 0.16 and a eee 0.60. 
Wels Long-wave emission from the earth's surfaces (Ry) 


The surface of the earth,:as any other body, radiates energy. 
Because of the earth's temperature most of the emission lies in the 
invisible infrared region of the spectrum. The earth is commonly 
assumed to emit and absorb energy as a gray body in the infrared region. 
The formula for the earth's surface emission may thus be written as 


- 4 en 
Rf sheavat T2273) (4.2) 


according to Stefan-Boltzman Law 


where o is the Stefan-Boltzmann constant 


= 5.670% 10°° Wonca ie ot 


qT, is the temperature of the surface in °C and ES is the 
infrared emissivity of the surface. The emissivity of a black body is 
unity. Typical infrared emissivities, expressed in per cent for the 


various surfaces are given in Table 4. The data are from Sellers (43). 
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TABLE 4. Infrared Emissivities (e,) 


(Per Cent) 


Surface 


fresh 


Snow, ice granules 89 
Ice 96 
Soil, frozen 93 - 94 
Sand, dry 84 - 90 
Sand, wet 95 
Gravel, coarse 91 = 92 Z 
Ground, moist, bare 95 = 98 
Ground, dry, plowed 90 

' Grass, high dry 90 
Field and shrubs 90 


Pine forest 


As indicated in Table 4 the values of E. for the various 


surfaces vary from 82 to 99.5%. 
4.2.4 Long-wave counter-radiation from the atmosphere Ry) 


Although the atmosphere is nearly transparent to shortwave 
wave radiation, it readily absorbs the terrestrial radiation. The 


principal absorbers being water vapor, atmospheric gases and clouds. 
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As a result only about 9 per cent of the terrestrial radiation escapes 
directly to space. The atmosphere in turn reradiates the absorbed 
terrestrial radiation, partly to space and partly back to the surface 
(counter radiation). The counter-radiation which is absorbed by the 


surface Ry? for clear skies is 
4 
= ate 2 
R ‘ ee 08, 273) (4.3) 


where E5 is the effective emissivity of the air and has been approximated 
by 


-—ce 
€ = a-— be a 
a 


a, b and c are empirical constants and oe is the vapor pressure 


in millibars. 


Geiger (45) uses 


"a = 0.82 
b = 0.25 
c = 0.094 


The vapor pressure a has further been approximatedby (Apprendix IV) 


_ wees" 4 7 2900/ (1,273) (4.4) 


where w is the relative humidity (Appendix V). 
The counter-radiation Ry from a cloudy sky is expressed in 


terms of the value for a clear sky (45) by the equation 


Se corr) = (14Mw ) e4 (4.5) 
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where M is a constant and We is the extent of cloud cover (Appendix 
VI), which for overcast skies is.1.0 and for clear skies is.0.0. 
Radiation from the clouds, as indicated in equation (4.5), does not 
increase linearly with Woe but almost as 2 quadratic function. Other 
than depending on the extent of the cloud cover, counter radiation also 
depends on the height of the cloud, as indicated by the values of M 

in Table 5. In general, the higher and thinner the clouds the less they 
will contribute to the counter radiation reaching the ground. Typical 


values for M are listed in Table 5 (45). 


TABLE 5. Values of M 


Cloud Type Height (m) M 
Cirrus 0.04 
, Cirrostratus 8,390 0.08 
Altocumulus | 3,660 O.47 
Altostratus 25440 O220 
Stratocumulus 15220 0.20 


Stratus 


46955 Sensible Heat Flux (S) and Latent Heat Flux (L) 


The temperature of the earth surface is usually different from 
the temperature of the lower layer of the atmosphere. Consequently, a 
vertical heat flux exists between the underlying surface and the atmosphere 
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by the turbulent heat exchange in the air layer near the ground. 

The determination of these vertical turbulent heat fluxes usually 
presents the greatest difficulties in any heat balance problem. Because 
of the complexities involved in developing suitable equations for these 
two fluxes,detaileddescriptions will not oe presented here, however they 
are available in Sellers (43), Budyko (44), Geiger (45) and Matveev (47). 
The most common method for obtaining the formula for the vertical 
turbulent heat fluxes in the air layer near the ground is to assume 

that the process of turbulent diffusion is similar to that of molecular 
diffusion. 


Hence for sensible heat flux S, the equation is 
S=- pe (AT/Az) (4.6) 
and for latent heat flux L, the equation is 
Boek res, 
Ti pL K (Aq/hz) (4.7) 


K and K are the corresponding eddy diffusivities and q, T and z 
are the specific humidity, temperature and height. Assuming S and L 
do not vary within the air layer thickness of one or two meters above 
the underlying surface, equations (4.6) and (4.7) can then be integrated 
from the surface to a height z, ome or two meters above the ground to 


give 


io) 
HT} 


pe D, (T.-T,) (4.8) 


and 


& 
i 


pLD (q5-4,) (4.9) 
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or 


L = (0.622pL.D_) (e,-e,)/P (4.10) 


Further, if it is assumed that the diffusion of heat takes place by 


the same mechanism as the diffusion of water vapor, De = Di» hence 


equations (4.8) and (4.10) can be written as 


S = A(T.-T,) (4.8) 
and 
L= B(e.-e.,) | (4.10) 
where 
A= D 
BS 
B= (0.622pL D)/P 


D; Gees ae P are the diffusitivity coefficient, surface 
vapor wtacees air vapor pressure and the total pressure. 

From Appendix V, the saturated surface vapor pressure e. is 
given by 


= 923-4 ee UE) (4.11) 


ss 


and that for the air vapor pressure is 


wero’ 9 2200/ (Tat273) (4.12) 


a 


For the top surface with a certain moisture content Ms> the surface 


vapor pressure e. is 
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Pes) 5900) (157272) 


e =M + (1-M.) e, 


II 
K< 
io) 
+ 
4a 
T 
<4 
o 
@ 
re) 


(4.13) 


D, the diffusivity coefficient for heat and moisture in-the air layer 
near the ground are presented in various forms by different authors 
(43, 44). The principal properties of the diffusivity coefficient 

are that it changes considerably with the wind velocity, increases as 
the surface roughness increases and that the coefficient is higher in 
dry region. It also will depend on the stability of the air layer. 


From Sellers (43), the diffusivity coefficient is given by 


Di=?f i ou (4.14) 


where D is inm./sec, g is a constant with values around 0.002, u is 
the wind veteeity in m./sec and the coefficient f is extremely variable 
with values ranging from 0.0045 m./sec to 0.0007 m./sec, depending 

on the type of surface. The coefficient -f increases as the surface 
roughness increases. 


Values used in the calculations are 


f = 0.003 m./sec 
g = 0.002 
“w= 4.5 m./sec 
426 Heat Flux to the Ground (G.) 


The insulating layer lying over the ground was assumed to 


have no heat capacity, therefore, with reference to Figure (29). 
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= G 
oe s-l 
rs Ty = qT, ar (d/k, ) Ce CaaES) 
or 
G, = (k,/d) (it) 
LOND ey | The insulating layer (U-factor) 


The layer was assigned a U-factor, k,/d, having an average 


winter and summer values, ue and Us. 


d/k, = d_/k Or42 71-1 =. -0..83 
ef Ss = = 
=d_/k 0.42 > tn) > 0.783 
w Ww 
“4.2.98 Energy Balance 


Net: heat flux to the ground surface is 
F = (1-r) Ro ap Ry - Ris - §=L (4.16) 


For energy balance at the surface, F is also equal to Go» the heat 


flux to the ground. 


4.3 The Input Variables Required for the Solution to Equation (42) are: 
(a) Radiation Parameters: 
Ros AR, Pps ae a; b, ¢ 


(b) Relative Humidity: 
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(c) 


(d) 


(e) 


(f) 


(g) 
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Cloud Cover: 
Wo M 
Sensible and Latent Heat Parameters: 


0, Lo» eet ae 5 ells o 


Air Temperature Parameters: 


Ti AE Oe 


Surface Conditions: 


ye da» di US to» Lor ki? af 
Ground Parameters: 
ke, Ke? PCE» PC if? Pkg eS 


Standard Input Data Set (Norman Wells) 


i) 


id) 


Radiation Parameters: 


Ro - 100 WebESy a 


AR. - 100 MALES IAC 


$2 - 70.223 
€ - 0.90 
s 
a - 0.82 
b — (eZ) 
c - 0.094 1/mb 
G == 5.67 "x 10° watts/ fe 0x4 


Relative Humidity: 


Ww - 0.65 


ie 


@ seta 


i 
= 


 . 
ac 


et ee ee 


a ae Ps) i J 
a. ¥ * ee bite, a / 


* Ue, & 


ea i Ane,0 “ is 


"Hawt + 8 


1 Ks “ - LX eae 
. 
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, = = i” 
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iii) Cloud Cover: 


iv) 


v) 


vi) 


Ww E> (0S 


Sensible and Latent Heat Parameters: 


ooh = T ete/n 

Le 22315 ‘Ne Joules/kg 

P = = LOLS smb 

f - 0.003 m/sec 

Cee e000 

u - 4.50 m/sec 

ae bedi 10° Joules/kg °C 


Air Temperature Parameters: 


(iy ea eee re 
a 
AE 22206..°C 
a 
nr £0.29 


Surface Conditions: 


M, - 0.30 (surface moisture content) 
r = (0.16 
s 
1s - 0.60 
Ww 
d = 70525. Ml: 
Ss 
d=) 0.38 m, 
Ww 
k. = 0.30 watts/m°C 
k. = 0.17 watts/m°C 
ty =e 10.42 
ef. *=8-0583 
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vii) Ground Parameters: 


Moo - 0.30 (soil moisture content) 
ke 1.70 4watts/m°c 
kY.— DL.10°watts/m°c 

uf 


pc, - 1 OO es 10° Sonles/me 26 
PC em 2.40 x ane Foulee ine. oc 


Deas, bs 2 OL 10° Poulenc 


4.5 Introduction of the Analytical Section 


For a given set of input variables, equation (4.16) along 
with the equations of heat transfer in the ground (Part III) can be 
solved for the ground temperature variation throughout the year. This 
will be done numerically in the following section. However, because 
of the large number of input variables involved and the uncertainties 
in some of the empirical expressions for the heat flux, the absolute 
determination of the ground or surface temperatures may be questionable. 
Qie can however predict the changes in the ground and the surface 
temperatures that would result from changes in any one of the controlling 
parameters. This approach which is essentially a 'sensitivity analysis' 
will be taken with the model presented. First this will be done 
analytically. The analytic calculation will require the neglect of 


the interaction of the ground and the surface heat fluxes. 
(ee roa Analytical Section 


The main objectives of this section are to obtain analytic 


estimates of the effects of the various heat fluxes acting on the 
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surface of the ground cover. These effects will be expressed in terms 
of an effective surface temperature, an effective "heat transfer 
coefficient", and an average influence coefficient. The values of these 


parameters will be compared to the values obtained numerically. 


Ae 5 2 Effective surface temperature and effective "heat transfer 
coefficient" 
A linear approximation to equation (4.16) can be obtained by 


expanding it in a Taylor series about some temperature eee Thus 


F(T.) = F(T,*) + 9F/aT, (iat yt (4.17) 


Tc, 
Ss 


where terms in a higher order in the temperature difference Ciset*) 
have been neglected. 
ig Tee is chosen such that F(T. *) = 0, then equation 


(4.17) becomes 


F(T.) = oF/9dT (Tee) CIS) 


| 
s| 
T * 
s 


The term R/O, | in equation (4.18) is essentially an effective 
IPs 
heat transfer coefficent which could be called Doce? thus giving 


a -T * 4.19) 
F(T.) hog 1, ee ( } 


The temperature T 3 is the temperature the surface would 
have if there was no heat flux to the ground. It will therefore be 
close to the true surface temperature when there is a thick insulating 


layer of snow or vegetation on the ground. More specifically one would 
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expect equation (4.19) to be a good approximation when he 2>1U3 


Ff 
where U is the U-factor or conductance of the surface insulating 
layer. This approximation will be checked in a latter section. A 
temperature similar to i when used in heat transfer literature 

on human comfort or other similar topics is often called an effective 
or radiation temperature. The former term will be used in this 
thesis. 


To find To equation (4.16) is set equal to zero and may 


be rewritten in the form 


4 -5900/Y 
e 


HS GN came (4.20) 


where Y = Tf. * + 273 and A.» Ay and A, are known constants. 


> 
i} 


4 
HEA x (RE (isn) + oe o(T_+273) 
: (corr) 


629 L973 
+ (0.622 pL.D/P) Moe, + pc D(T,+273)} 


bo 
i} 


1/pe D x (e. o) 


23.4, 


> 
i] 


1/pc D x { (0.622pL_D/P) Me 


Defining the functions 


: 4, .-5900/¥ 
Pope ee A Ae 


3 -5900/Y 


2 
£"(Y) = 1.0 + 4A,Y + (5900 A,/Y") e (4.22) 


ok AS JW dee as Et @altdl 209 
wy 2 Nom, i 


‘ 7 A 
Ate Leer aa. IS GO. OE 


SP ale maei ‘a a rAih4 ad % 


ys 4 , SHO AT oni Ty 


oy Psa 1 
Ee ee 


io hed hte ney i 


‘, Tyaparts-- oe La iw f 


4 
PP owe § 


7 i . oe sa 
+ oh LP Rca ey * iP 

7 rh) ie Lore ” : - 
ae he, sn-78 om. xe F 4 Berean ay ern { hr ~ + * 


os ; 
AS D0 Oy » 
oe , aL 7 iy ; ’ 
. M, 7 -? a“g J i vps Ly 4 3! ~ a. mye A ‘b 
r Ls : 7 hy F ) 
bug onm ak TOC ggg aE Ye 
Lge os ol (20 NY Rae RR RAR 
i Fee Oy ee : J aa ae 
j Rs -s me i ai a5 
= \ fav) a ay 
x} ‘ « . ie a : zw 9 = ‘5 ot 
: ; x ‘ zy ! ‘ A 


if 
are 


< ’ i € 4 1) i 7 i & Li a 
Lee ay Lene ae th eu ce! phabs r* Mids te ay ™ 


fee 


P 
' > 
\ 
U uv 
+ 
; 
bil 
—— Ae 
- 
ot 
1 
25 
* “a. * 
i 


BPO Pah RE RN 
' ‘ ‘ a4 eS b} ) 


ise oS - iw) pe tee ae we y heer 
ey ee ee ee ee 
Cc AGA Bg tat! 1 SM. a4 Pan 
yeas » " @ a i ee * - 
& ® "7. i ; i ic&h 
Saeteey = ; Pllete , Aye 
REE pd | Cee Ny 
: . (4.50) 2 U agit 
ee ( - Cn yee om ny y My i 
Aoae coal . as F rn & Y | 7 ; 
i = ve 4 Sto) } e sq i ~ yA 
- r a | . 
a : ir oa - tre 
Piira> 7 |F ea’ ia anbakionl 
E; ; , is * 
oan: 
S ei Te, ap bs 4 
" q rs Poe ae i 6) er aS ee cy) q 
a MI eae : ain 
‘ane? Ne ; ‘ : : - } 
4) wee pe : J 4 Vp Aa ¢ a €y¥) 2 ; 
= ee | : =n 
was i 7 s ‘ 
" ‘ 1 7 % 
a 4, fe { : ; i 


_T E a, 


s 


98. 


Equation (4.20) can be solved by first assuming an initial 
approximation to a root, then generating successive approximations from 


the iteration 


ze fe ' ’ aes 
an Y, £(Y,)/f (Y,) (Newton's Method) 
ge for the entire year for the standard data, given in section (4.4) 


was then obtained as tabulated in Table 6 along with the values of air 


temperature, T.: 


TABLE 6. Values of i for the entire Year 
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It will be noted that tee is lower than qT in ae and 
higher than qT. in summer. This difference may be attributed to the 
fact that the value of the reflectivity of summer vegetation is very 
much lower than the reflectivity of winter snow cover, the respective 
values being 0.16 and 0.6. The greatest difference between tee and 
T occurs in the winter when the additional heat flux due to solar 
radiation is minimum. 

To find the effective "heat transfer coefficient", equation 


(4.16) must be differentiated with respect to Tine This gives 


3 
OF ie 9 - 
oF/3T. 4 E. o(T,+273) eee 


- (0.622 pL_D/P) Moe" (5900/ (1_+273)7) 972 900/ (T.,+273) 
Hs . 
hore aH,/ 27, | (4.23) 
T * 
s 
The effective "heat transfer coefficient" h is thus the sum of 


eff 


the terms due to radiation, sensible heat and latent heat. These terms 
are labelled hes he and hy respectively and values are given for 
various times during the year in Table 7. It will be noted that hy 


and h, are time dependent while he is not. The effective "heat 


L 
transfer coefficient" Doge that results, ranges from a low of 18.48 
2 F 
ate in winter to a high of 30.0 watts/m-°K in summer. Its 


; 2 
average value for the entire year is 22.51 watts/m-°K. Values of Hee 


obtained from William (51) for latitude 40°N to 60°N indicates that 
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for average conditions of exposure, the coefficient varies between 
19.39 to 29.09 meceein ok (16.68 to~25.02 keal/m’h °K). This analytical 
value of effective "heat transfer coefficient" Doge? will be later com- 


pared to the value of he obtained numerically. 


EE 


TABLEL?7. Values of hs for the Entire: Year 


ses os “yy eee 


(L-w radiation) (Sensible) (Latent ) nacre ook 


45. 3 Analytical average influence coefficient for the standard 
data set 
The main objective here is to find the change in T, caused 
by a change in x; where Xs is one of the controlling input parameters. 
If the net heat flux to the ground is zero, then equation (4.16) 


becomes 
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if shee 4 4 
0 = RQ t) E. o(T.+273) + E.€, o(T +273) 
(corr) 


~ es a Y, - (0.622 pL D/P) (e,-e,) 
which implies that 
£(T., x,) = 0 | (4.24) 


Taking the total derivative of equation (4.24) with respect to the 


independent parameters Xs gives 
= = £ | e 
df /dx, 0 of/9T. aT /dx, + of / 3x, 


implying 


dT. /dx, oo Ne {df /dx,/3£/3T.} (4.25) 


of/oT. was obtained previously in finding the effective heat transfer 


coefficient. The average influence coefficient is 
dT / dx, 


which is equalled to 
1 . 
f dt fdx. dt (4.26) 
Se aL: 
0 
In integrating equation (4.26), the nominal values of T, is assumed 


to be TS*; the value which gives F= 0. Taking the partial derivative 


of £ with respect to the independent variable x» of /dx,: 


2») d£/ dR, = (l-r) 
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dR /oR = 1 
Ss s 


OR/ 9 


(AR/R_) = R, sin(wtt,) 


. a£/9R_ = (1-r) 


9£/3(AR/R.) = (1-r) R. sin(wt+$,,) 


44); O£/ 0x = ane 
iii) df/de = -o(T 4273)" +e fone 4a b 4273)4 

s s “a a 

corr 

iv) 0of/du = mite eae sb ~ (0.622pL_8/P) (e,-e.) 
v) of/ ow = Ze OnMews ce meek 4273)" 

c s cya a 
vi) 9£/ 0M, = - (0.622pL D/P) (e,.-e.) 
vii) df/dw = es OCT 4273)" (1+Mw *) (bee e 4) 

8 oes c 7sa 
aF (0.6220L,D/P) M, ey 

viii) d£/9T = Gee CE 4.273) >+ e o(14+Mw 2) (7 4273)" x 
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aT Lot a sul: 
a a 
oT, /o(AT/T) = qT, sin (t+$,,) 


- of/9T = Of/9T x dT /T 
a a a a 


and 


o£/3(AT/T) = 3f/9T x oT /9(AT/T) 
: a ane ie 4 
ix) of/ oT. = =e o (T,+273) ee 0.622oL D/P X 
2 
M e _5900/(T_ +273) 
s'-ss Ss 


As stated earlier dT /dx, can be determined by using the 


chain rule 


dT, /dx, =- {8f/ dx, /3f/oT 5} 
The average 
al 
dade) Viel aig dx ads 


where it is assumed that its = z,*- 

The integration of dT, /dx, over the entire year was done 
by trapezoidal rule. 

Define the influence coefficient to be the change in 
temperature D OL T, caused by a percentage change in the input 
parameters. Note, parameters fr, Eg Ws M, and w which represent ratios, 


will be expressed as percentage, 1 to 100, for the purpose of uniformity 


in the defination of the influence coefficients. That is for instance, 
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the standard E. used is 0.9 which will be 90% for use in calculating 
the influence coefficients. The one exception is the air temperature. 
Its influence coefficient will be expressed in terms of a temperature 
change in °C of a or es for a given temperature change in °C of 
it 

The influence coefficients will be used later to determine 
which are the most important parameters in so far as their effects on 
the ground temperature are concerned and also they will be used to 
estimate changes in ground temperature that would result from given 
changes in the input parameters. Ideally, to be of practical use, the 
influence coefficient for a particular parameter should be a constant 
over the range of that parameter and be unaffected by changes in other 
parameters. That is the effects of a parameter should be linear and 
uncoupled. As'can be seen from the equations in this section (4.5.3), 
for the various influence coefficients, this cannot be exactly true for 
any variations of the input parameter; however, it may still be a 
reasonable approximation for small changes in the parameters. 

To get some estimate bf the likely values of the influence 
coefficients and the range of their values, Table 8 gives the average 
coefficients for the two data sets. One is the standard data set from 
section (4.4) and the other is a modified data set obtained by changing 
selected parameters within the limit of their reasonable ranges. It 
can be seen that most of the influence coefficients are the same for 
both data sets, however, some do vary. In particular the coefficient 


for (Au)/u changed very significantly between the two data sets. This 


is because it depends on the temperature difference, Cus wer which in 
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turn is affected by all the other parameters. The effect of the average 
wind velocity u, is thus strongly coupled. The influence coefficient 
for M, changed by about 30% between the data set. This could be 
attributed primarily to strongnon-linear dependence of vapor pressure 


on temperature. 


TABLE 8. Analytical Values of the Average Influence Coefficients 


ee ee, 


Input Average Influence Coefficients 


(4) 


Parameters Standard Data Set Modified Standard Data Set 


T 


iS} a 

0.024 0.025 

ian ie 0.090 0.070 
A(r) -0.051 -0.045 
Ae.) "19.047 ~0.047 
(Au) /u ~-0.028 -0.043 
ime)? 0.018 0.017 
Aw.) 0.027 0.026 
act.) ° -0.017 -0.023 
i {0.819} {0.812} 
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(2) R,=w 
(3) M. = XMS 
(4) some of the parameters of the standard data set are modified 


within their limits from their individual nominal values. 
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4.6 Numerical Model 


Basically the numerical model to be presented is similar to 
the numerical model presented in Part IJ of this thesis, except for the 
different boundary conditions and changes in the definition of the 
various dimensionless variables. The problem now has, instead of an 
imposed temperature Sates condition, a varying heat flux at one 
boundary. The changes in the definition of the various dimensionless 
variables are necessary so as to accommodate the new parameters. 

The varying heat fluxes are those discussed in the early 
sections of this part (Part IV). The insulating layer is assumed to 
have negligible heat capacities and values of the insulating values 


are from Appendix III. The parameters in the ground remain unchanged. 
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Governing Equations 


CG. = (l-r) Ro + Ry = R 4 So = lb (aa 27) 
sa sh 7 G. 
ig = qT, ay d,Go/k, 
dx dh) /at = CG. | - GC, 
Ax oh, /dt = Cla = Cian 
h = pc-T i< qT. 
h = pc el ap pL, laa T, 
Ge ee Ca 
where k = ke 1 es Ty 
= k if le= Ty 
AeOs2 Normalization 


Defining the various dimensionless variables. 


Writing 
ar Re /ou, where eae 1 year 
AT «= RLx,/k¢ 
6 = (1-T,,)/AT 
¢ = h/pl, 
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(4.28) 
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9 iy = R ' ~~ — ~ 
(x AtpL,e/t,) 86,/3t = R, G)_| - KAT(@,-6,)/x Az 


Ag, = (on At - k'(6,-8,)} Atx2/0" i=1 
Ag, = {k"(0,_4-8,) ~ K'(,-8,,,)) Ar/ac? av 
o* = 5, 8 <0 $21 23 30 
ORS 80+ 1 o> 0 tai, 2s 5a% 
ki = k 6 41 18,70 
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Rewriting the governing equations in the normalized form 


Oh ae = ° asay ; 
R, 1 + AR/R. sin 2m (t+$,) (4.29) 
toe ee o(T +273) */R | (4.30) 
oh ¥ sa a Ss 
corr 
Py =iermo CT +273) "7R (4.31) 
cai + s S s 
i yee ee De c 
S A(T. TD/R, (4332) 
' = ! — /R ks 
L Ble. e/R, (4.33) 
tore _ ' We re es CoM ee Be - 34 
Ge = (1-r) Ro 7 RQ. \ Ro , S L (4.34) 
uy R 7 ' 3) 
oF 8) ot (ae ately c (4 ) 
nit ' ae te , i= 36) 
Ad, {G. Az - (8, 8,)} Atx2/Az i=1 (4.36) 
2 
Hie e 2 : ; 
m+1 m ™ 4.38 
ih eS | G38) 
6, =o, /S, >. 
(4.39) 
m+1 m+1 7 
a = (>, -1)/S ¢ >, 0 
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406.3 Method of Solution 


(a) Initially, set all the ground temperatures equal to an approximate 
initial temperature, 6J (ie. o = OS) )r. 


The initial enthalpies are then calculated. 


(b) Using equations (4.29) and (4.30) to compute R.! and ae 


respectively. 
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o(T 4273) °/R (4.30) 
s Mere a s 
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(c) To solve for = 


In equation (4.34), only R'4> Sandal are functions of To: 
On substituting these expressions into equation (4.34), the 


resulting expression is 
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Then substitute ce into equation (4.35) 


pol 2D ’ \ 
oe ct + (R.d/k, AT) Gs (435) 


and writing ae = (T-T,)/4T, equation (4.35) becomes 


ere rae te (4.40) 
where Y = Ie +8273. <and A> Ay» A, are known constants. 
A, = {(1-r) R. + Ese, o(T,+273)" “+ penal) 
corr 

+ (AT (8, + 273))/d/k, + BMe. /(oc Dt1/d/k, ) 
A, = €,o/ (oe Dtl/d/k; ) 
A, i Bye" "/ (oc DEL/a/k,) 
f(x) =¥-A + AY +A, ao (4.41) 
Deas gee 4A,Y ‘ 590A, /Y" Pee: (4.42) 


Equation (4.40) is solved by first approximating a root, then 


by means of successive approximations from the interation 
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With the above T,» values for Rye S' and. L" are obtained 


from equations (4.31), (4.32) and (4.33) respectively, 


foe 
! = 
Rees 0 (T,+273) /R. (4531) 
od = re os 8 = f 
Ss A(T. T DIR, (4.32) 
hin a Dp 
i, B(e. e D/R, (4.33) 


Then compute G.! using equation (4.34) 

' = <= ' LF oD es Ls pe sei (al Mie ee a pel | 
G, (1-r) R, ze R. : R. 4 S L (4.34) 
Compute the enthalpy changes at each nodal point, using equation 


(4.36) for the top layer and equation (4.37) for all the remain-~ 


ing ‘layers. 


agy = {C,'A¢ - (6,-9,)} Atx2/ Ace i=l (4.36) 


£(8, 79,) + (i4,794)) oie eR (4.37) 
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Then using enthalpies from the last time step and plus these 
enthalpy changes, new values of enthalpies at each nodal point 


are calculated by using equation (4.38). 
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(g) The enthalpy o is then used to calculated pats by 
using equation (4.39) 
mi] m+1 
ew ee ae” 
(4.39) 
mt+1 
~ AG Nesey sae $, > 0 
aes 1 
e and ae are the new enthalpies and temperatures one 
: m cage ; 
time step ahead of OF and e. respectively. 
(h) Proceed as before, as the time step advances, the initial 
a 
enthalpies and temperatures are those just computer, oo 
m+] : = : ‘ 
and | respectively. The computation is over when the 


temperature distribution for an entire year approaches a steady 


state. 


4.7 Results of the Numerical Calculations 


As stated in Part II of this thesis only the steady periodic 
solution will be of interest, the transient behavior produced by starting 
from a constant temperature will not be considered. Usually it was 
necessary to carry out the caiculations through only five yearly cycles 


of temperature variation in order to obtain the steady periodic condition 


(refer iiosrart 21). 


[Ned / elk Heat Fluxes 


Figure (30) shows the values of the various Heat Fluxes and 
also the Net Heat Flux to the ground surface during the fourth to the 


fifth year. As indicated, all the heat fluxes peak at about the same 
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FIGURE 30 — Values of the various heat fluxes to the ground 
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time of the year, that is, when the air and surface temperatures are 
maximum and correspondingly as the air and surface temperatures are 
minimum. The net incoming Solar Radiation (1l-r) Ro» which is not 
directly temperature dependent, however peaks during the summer months 
and has a minimum value during the winter. All the component heat 
fluxes have positive values, that is they are in the direction defined 
as Figure (29), throughout the entire year except for the Sensible Heat 
Flux, which varies from a minimum negative value to a maximum positive 
value. This is because during the winter, the surface temperature ie 
is lower than the air temperature qT. and in summer, Ty is greater 


than To The empirical equation for S is 
S = pe D(T-T.), 


as such §S is negative in winter and positive in summer. 
The net heat flux to the ground on also varies from a minimum 


of ~12 aeee ne to a maximum of 15 watts/m*. The empirical expression 


Ch Cae em ae esta Say! 


The negative values of GC. during the winter months indicate that heat 
is given off by the ground during the winter, and similarly positive 


values of G indicate that heat is absorbed. 
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TABLE 9. Comparing the Various Heat Fluxes from Various 


Sources for Norman Wells (65°N) 


Data Source Net Radiation |} Latent Heat! Sensible Heat 


L S 


2 
watts/m Waeea ya 


Sellers (43) 


Mean values for 60-70°N Tae 
22 Mare, (52) _ 
3. Standard Data Set -22.4 


4. Modified Standard Data Se 


As indicated in Table 9 the vaiues of the various heat fiuxes depends 
very much on the values of the input parameters. The valves from 
Modified Standard Data Set.agree favourably to those obtained from 
Sellers, whereas those values from the Standard Set are significantly 
different. The modified set are just the values of the standard set 
with some variations in some of the controlling parameters such as 
TR 5 AR. Ee» Spring .a Balt, d. and di. The variations are however 


a S 


within the range of variations of each of the individual parameter. 
Quien Temperatures Distributions 


Figure (31) shows the values of the various temperatures, 
namely, us (air temperature), Te (surface temperature), Te (ground 
surface temperature) and Sas (temperature just below the active layer) 


during the fourth to fifth year. The graph shows that in winter the 
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FIGURE 31 — Profiles of temperature distributions with varying heat fiuxes at the 
surface. U = 0.76 watts/m2 — °K, (U, Uy )/U = 1.0 and ¢,, = 0.0 
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temperature deep in the ground ae is warmer than the ground surface 


temperature cod which in turn is warmer than the surface temperature 


T.: The air temperature qT, is cooler than the ground surface temperature 
Oe but is warmer than surface temperature TS: During the summer months 


the reverse is true, that is To Lsncooler than T which in turn is 


gl 
cooler than T 5° The air temperature qT, is warmer than at but is 
cooler than TS: It also will be noted that the averages of these 
temperatures over the year are not the same. The differences in the 


average temperatures between the various levels is due to non-linearities 


in the equations. 
x = 001° 
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The difference in ( ) was studied in Part II and the difference 
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section because of the non-linear nature of the heat flux equations. As 


indicated in Figure (32), the average temperatures are 1 = -6.22°C, 
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i 6ige, Ty = 73+93°C and ae = ~4.33°C. Note that (T.-T),) is 
by far the largest of these differences indicating the importance of 

the insulating layer in the determination of the net difference between 
the ground temperature and the air temperature. Measured values obtained 
by Brown (29) for the difference between the average air temperature ts 
and the average ground temperature ae (50 to 100 feet deep) varies 
trom 2.7 to%4.4°CG (or -5 to 8°F). However in one of his other papers 
Brown (50) obtained another set of values of ve - To (12" deep) for 
Norman Wells which ranges from 0.7 to 11°C. The measurements were made 
under different conditions at different sites. The computed rote for 
CE Stee)” eS 1.9°C. These large variations in the value of (Ton ) 
are mainly due to large number of varying influential parameters which 
the vaiue of Cae) are dependent on. Factors that are particularly 


influential are net radiation, vegetation, snow cover and ground thermal 


properties. 
NT he) Depth of Thaw Penetration 

The depth of thaw penetration in Norman Wells is about 0.625 m 
(approximately 2 feet), as indicated in Figure (33). The result compares 
favorably with those obtained by Gold (33). 
4.7.4 Average Influence Coefficients 


Earlier in section (4.5.3) the analytical model for finding 
the average influence coefficients was developed by assuming that there 
was no interaction between the ground and the surface heat fiuxes. The 


values for the standard data set were presented in Table 8. 
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In this section the values of the average influence coefficients 
are obtained numerically and the interactions between the ground and the 
surface fluxes are taken into account. The values are obtained by 
varying a single parameter at a time and then divide the change in the 
average surface temperatures Ge and the change in the average ground 
temperatures (T we) by the change in the varying parameter. The range 
of variation for each individual parameter is a carefully chosen range 
within the realistic climatic conditions of the location under study. 

The values of the average influence coefficient for the standard 
data set are presented in Table 10. The table also shows that the set 
of values obtained analytically compares favourably with those obtained 
numerically. From these coefficients and an estimate of the possible 
range of each parameter (column 5) the magnitude of its effect on the 
ground temperature can be estimated (column 6). The crder of importance 
of a parameter in determining the difference between ave and cM or 
art) for a given ie will be indicated by the magnitude of OT os 
in column 6 of Table (10). That is to say the parameter inducing the 
largest a: will be the most important and so on in descending order. 
The order of importance are (AU) /U, (AU_)/U_; Moo? (AV)/V and 
(AR, )/R, and is in rough agreement with the field observations summarized 
by Brown (53). 

The range of ee possible for one air temperature is 16.77°C. 
(sum of ae in column 6). The range of (1-7, ,,) (12 inches deep) for 


Norman Wells obtained by Brown (50) is from 0.7 to 11°C. This can explain 


the occurrence of isolated pockets of permafrost in regions as far south 
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TABLE 10. Numerical Values of the Average Influence Coefficients 


Average Influence Coefficients Range of 


(Analytical 


Input 


+3 | 
4 | 


Parameters 


(AR_)/R. 


A(r,) 0.66 
A(x) 0.30 
A(e.) 0.40 
A(R) 0.99 
*A(M,) 0.48 
A(w.) 0.42 
(AV) /V 1.26 
«(AU ,/U, 3.60 
* (AU) /U,, 5.60 
{AT} 8427 
Neues hae’ 0.60 


*act,) 7? 


- density of soil 


soil 


Qa M - soil moisture 
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as Edmonton and isolated pockets of thawed soil as far north as Fort 
McPherson (67°N). 

Most influence coefficients were found to be constant within 
20% over the range of the input parameters. Typical results for the 
non-linear cases are shown in Figure (34) and Figure (35). The applica- 
tion of these influence coefficients will be discussed later in section 


es ek) is 
Genks D Effective ‘heat transfer coefficient' 


In section 4.5.2 the analytical value for the effective ‘heat 
transfer coefficient’ was determined by differentiating equation (4.42) 


with respect to the Ty to give equation (4.23). 


hoce = oF/ dT. (4.23) 
4 T * 


s 
The value of Doct was then determined by assuming the nominal value 
of 1 to be T*> the value which gives F = 0. 
In the numerical model, the effective 'heat transfer coefficient' 
h is determined by dividing the change in the heat flux to the ground 


eff 


G. by the corresponding change in the average surface temperature T.: 


These changes can be brought about by gradually increasing the heat flux 


from deep in the ground. 
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SURFACE MOISTURE CONTENT, m,,% 
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AVERAGE GROUND TEMPERATURE ( 
AVERAGE SURFACE TEMPERATURE 


FIGURE 34 — Variation of M, with T, and T,.. 
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FIGURE 35 — Variation of M,, with T, and Tg. 
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HEAT FLUX TO THE GROUND, G,, watts/m* 
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FIGURE 36 — Relationship between the heat fluxes to the es and 
the average surface temperature 
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Empirically the model can be written as 


-T *) (4.44) 


be = AG /AT 
Ss Ss 


From Figure (36) the value of he obtained by this numerical method 


ff 
is 24.24 Sear which is slightly higher than the value of Doge 

(22. OL We ete ae) obtained analytically, However both the values are 
well withinythe range of 19.39 to 29.09 See mes ke obtained by Williams 
(51) for latitude 40° to 60°N. This also indicates. that the analytical 
method used is good for a rough approximation, provide the nominal value 


ale 


of T_ vis nearly equal to Te 
LEAL) Coupling between Air and Ground Temperatures 


Figure (37) shows the various temperatures at the various time 
intervals during the fourth year. Unlike the case in Figure (31) where 
the average U-factor U was 0.76, the U in this case is relatively 
large, about 6.0. Comparing the two graphs show that except for the 
ground surface temperature ae all the other temperatures namely, 
air temperature tT. surface temperature T; and ground BEIGUS isthe 
T...» exhibits the same behaviors. Pi. for this case closely follows 


ee whereas in Figure (31), ve closely follows ibs This means that 


for U = 0.76 the major thermal resistance between the air and the ground 


is the insulating cover and for U = 6.0 the major resistance is provided 
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by the ground itself. The average surface temperature T. remains 
unaffected by the change in U, however the average ground temperatures 
above and below the active layer, ae and es become colder, from 
-3.9°C and -4.33°C to -7.16°C and -8.31°C respectively. 

Figure (38) shows more clearly the effect U has on the 
average temperatures of the surface and the ground. For small U-factor, 
thus high thermal resistance between the surface and the ground, the 
two ground temperatures real and Bh are almost the same. This occurs 


because the amplitude of variation in temperature eal is small and 


thus no difference between at and ane can be created by the mechanisms 
discussed in Part II. In this case the ground temperature will be 
determined by the characteristics of the surface cover and the climate 
parameters and will be unaffected by the soil type of the moisture 
content. WisomoLe that the average ground temperature is higher than 

the average air temperature in this situation. This is because of the 
differential heat valve effect of the insulating layer, as discussed in 
Part III. As U increases, that is as the insulating properties of 

the ground cover become poorer, the ground temperatures decrease and the 
difference between Ty 


phenomenon. First, as U increases, the effectiveness of the insulating 


and T increases. Two factors cause this 
eke 


layer as a differential heat valve decreases. Thus the temperature dif- 


ference (T =f decreases resulting in the lowering of the ground 
Ss 


gl? 


temperature. Secondly, as the U-factor increases, the amplitude of the 


temperature variation in CES will increase. This increases the 


effectiveness of the heat valve effect caused by freezing and thawing 


of the soil and generates the difference between Ty and a 
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FIGURE 38 — Variation of U with T 
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An average U-factor of six corresponds roughly to a surface 
with compacted snow, 10 cm deep in winter and a 5 cm of vegetation in 
the summer (Appendix III). For surfaces with less insulating value 
than this, the effects of the insulating cover on the average ground 
temperature will become negligible. In that case the ground temperature 
will be controlled by the soil properties and the climatic conditions. 
Figure (39) shows the relationship between (U,-U_,)/U and the 
various average temperatures for a given U. The figure shows that the 
average surface temperature ap is independent of (U,-U_)/U but the average 
ground temperatures ae 


gl 


increases. The difference between Te and ae TSG Ce hor. (U,-U_)/U = ioe 


and any increase very significantly as (U_-U_)/U 


This indicates that variations in the U-factor of the surface cover 
are the major cause of difference between the average air and ground 
temperatures. | 

It will be noted from Figure (38) that the average surface 
temperature is relatively insensitive to the U-factor of the insulating 
layer. This implies that the surface temperature is controlled by the 
balance of the heat fluxes that occur to and from the air and that the 
heat flux to and from the ground ieeeos small to have a major effect. 
This is basically the assumption used in deriving the analytic approxima- 
tion in section (4.5.2). As stated in that section the requirement for 
application of that assumption is that U << hore: The value of U 
used is 6.0 watts/ Teer and the calculated value for Doce 1s) 24.24 


wates jee Listed in Table 1l are the values of qT, and Ts for 


an entire year. 
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FIGURE 29 — Variation of (Us — Uw)/U with T 
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TABLE LE Values on T, and ay for the Standard Data Set 


ah ft 


Ss s 
(Year) From Section U = 0776 U = 6.0 
(5 eo) (U,-U_)/U=1.0 (U,-U_)/U=1.0 
0.000 Seo =O =3.0.45.1 
0.125 -28.22 os Oo -27.88 
0.250 sol a fei, cy US TES =14.01 
C5375 See! 3.50 Big lal 
0.500 16697 16.42 UG da. 
0.625 1S bye: T3023 die Pals) 
0.750 - 2.03 = 1293 - 1.50 
02875 = Ake oD =20..0) =19.96 


1.000 ~31.37 ~30.81 | ESONST 


Comparing i to ites in Tabie uu shows that the difference 
between the two is very small, even when U = 6.0. It also indicates 
that for problems with small heat fluxes to the ground, the problem 
of determining the surface balance pennies uncoupled from the probiem 
of heat flow in the ground. Under these conditions, it may also be 
concluded that it is adequate to assume qT, = Tt for most cases, how- 
ever a better estimate would be on = T *. This is because by assuming 
Fi-= oo the effects of the heat fluxes are taken into account. 


Ss 


Even if U for surface cover is greater than Deep? T, may 


still be close to ube: because of the insulating value of the ground 


itself. A general test of the applicability of the approximation 
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T = ie is given by the equation, 
aa aie es 
i 3 Ae, 


where q is the instantaneous value of the heat flux to the ground. 
If the maximum value of the heat flux q is less than 24 eee iar 


the maximum difference of (Liste) Witl.~be Jess than 1°C. 
ay eal Applications of the Influence Coefficient 


As indicated in Table 10 the influence coefficients are 
relatively constant over the ranges indicated. Because they are 
uncoupled, the eienee coefficients can be used to estimate the 
effects of the disturbances which simultaneously effect more than one 
parameter such as disturbances to the surface cover. 

As an illustration of the application of the influence 
coefficients, the effects will be calculated for a disturbance to the 
surface cover which causes the following changes in the parameters: 

r_ changes from 16% to 8% 

e changes from 90% to 97% 

U_ changes from 1.2 to 2.4 (Gabieme=ek) 

and M, changes from 30% to 80% 

Keeping all other parameters constant, the values predicted for the 
changes in ay) and cee using the influence coefficients will be 


compared to the values obtained numerically. The changes in ts and 


Teo caused by the change in Ue and M, will be obtained from the 


graphs indicated as their relationship are nonlinear. 
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M, changes from 30% to 80% 


From Figure (34) 
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will cause a change in U and (U,-U_)/U. The change 


is tome. 7610. L.25. 


igure (38) 
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TABLE 13. Values of Computed T_ and Bias 


Original New Values (°C) 
Values (°C) Numericaily Computed by Using 
ie Mee Influence Coeff. 
T, wet ah -8.4 -8.4 
a -4.4 -3.5 ~3.2 


As indicated the values obtained for Bb. and et by using the influence 


coefficients and graphs compare favourably with those actually computed. 
4.8 Conclusions 


The problem of periodic heat flow to the ground in a freeze- 
‘thaw is very complicated because it involves about 39 to 40 parameters. 
The values of fhe heat fluxes to and from as ground are very sensitive 
as demonstrated by their significant difference between the standard 
data set and the modified data set in Table 9. Thus the absolute 
determination of the ground temperature is open to doubt, however it 
is possible to predict the change in the ground temperature that would 
result from a change in any one of the controlling parameters. It is 
however impractical to solve the whole problem everytime there is a 
change in one of its parameters. A much more practical method is to 
make use of the influence coefficients presented in Table 10. From 
these coefficients and an estimate of the possible range of each 


parameter, the magnitude of its effect on the ground temperature can 
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be estimated. The most important parameter is the parameter with the 
biggest effect. The order of importance of the parameters agrees with 
those obtained by Brown (53) from field observations. The most important 
parameter being the winter insulating layer, then follow the summer 
insulating layer. soil moisture content, wind velocity and so on. For 
instance the predicted quantitative effect of winter insulating layer 

on the ground temperature is +5.6°C whereas the effects caused by the 
winter reflectivity, r is 220587G. 

Two sets of values of the influence coefficients were obtained, 
one was obtained by the analytical method where the heat flux to the 
ground was neglected while the other was obtained numerically. The 
values obtained analytically for most cases agree with those obtained 
numerically. These values of influence coefficients, as demonstrated 
in section Roe could be used to predict the changes in om and ae 
caused by simultaneous variations in some parameters. 

The omni suggested that various approximations to the surface 
boundary condition could be applied. First and the simplest would be 
to assume qT = To This is a common assumption in some calculations 
of thermal regime in the ground. foe the standard data set used, it 
was found that the maximum value of |T, - 2 was 3.4°C. The second 
approximation which would include the effects of radiation, sensible 
and latent heat fluxes would be to assume se = peace fiat is the value 
that Ts would be if heat flow to the ground is neglected. This 
approximation uncouples the problem of determining the heat fluxes at 


the surface from that of heat flow in the ground. The sufficient 
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criteria for using this approximation is that U-factor << Deore: The 
maximum value of |T. - De was less than 1°C for U < 6.0 watts/m--°K. 
The third approximation would be the assumption that instead of an 
imposed surface temperature, a Neumann boundary condition exists at 
the surface. That is that the surface temperature is related to the 


heat flow to the ground by the equation 


as -T x 
q hoe MT, _ ) 


where both he and ia are time dependent. Both the analytical 


Bak 


and the numerical value of the effective ‘heat transfer coefficient", 


h obtained agree favourably with the value obtained from field 


eff 
‘observations by Williams (51). In this approximation, the very 
complicated heat flux relationships at the ground surface have been 


simplified by,a linear dependence. This last approximation would 


appear to have very general applicability. 
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APPENDIX I 


Data used for the Air Temperatures. 
Edmonton: 
From Figure (40), the average air temperature and the amplitude 


of variation for Edmonton are: 


Pasion 6 


ber 
i} 


AT = 15.96°C 


The values for Edmonton are obtained from the Annual Meteological Summary 
for Edmonton (31). Adopting the same procedure, similar values were also 


obtained for Norman Wells and Inuvik and are those tabulated in Table 14. 


TABIGH LG ies and AT for Norman Wells and Inuvik 


Amplitude of 


Locations Average Temp. 


T (°C) Variation AT(°C) 


Norman Wells 


Inuvik 


The values for Norman Wells and Inuvik are from Burns (30). 
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APPENDIX II 


Soil Properties from 'The Thermal Properties of Permafrost' 


by Virgil J. Lunardini (28). 


Fine Grained Soils 


TABLE 15. Values of ky and Ce for Frozen Soil, 'T = -3.89°C 


Moisture Content, Moo Thermal Conductivity, ke Heat Capacity, Ce 


(% of dry wt.) (watts/m-°K) (Joules /m>-°K) 


0.82 


TABLE 16. Values of KL 


and ¢ (2for Untrozen Soil, T = 4.44°C 
f£ uf 


Moisture Content, M 


he Thermal Conductivity, k 


Heat Capacity, c 
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Caiof ary iwti. ) 


(watts/m-°R) (ieules/m 43k) 
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TABLE 17. Values of Le 


Moisture Content, Mo Heat of Fusion, Le 


(Zoot dry wet.) Conies/m 


60 x 10° 
120 x 10° 


200 x neg 


APPENDIX: £11 


U-Factor for snow cover 


The average density of the snow cover is 250 als and its 


average thermal conductivity is 0.13 watts/m °C (Gold, 1963). 


TABLE 18. 


Year 


0.0 


0.083 


0.167 


0.250 


05332 


0.417 


0.500 


0.583 


0.667 


0.750 


0.833 


Osos) 


1..00 


Snow Depth (m.) 


Minimum Maximum 


Yon 
Depth Depth Uoy Ctin) 

226 0.46 
28 0.46 
«46 0.23 
-48 0.27 
-0 eo 
.0 co 
“(0 oo 
.0 G0 
0 oo 
.0 co 
.0 © 
3.5 0.87 
.28 0.46 


Snow Cover Depth for Norman Wells from Potter (37) 
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U-factor for Vegetation (Peat) Layer 


TABLE 19. Values of Kk and k 


£ £ for! eeat 


Dry Weight Moisture Content Thermal Conductivity watts/m°C 


eae ZL Of dry we. 


Values of the thermal conductivity were obtained from 


Lunardini (28). 
Consider only the conditions that give the largest and the 


smallest values of kf and ke for peat. 
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TABLE 20. Values of ue 


£ and U, 


nce 


Dry Weight Moisture Content Depth of Peat U-factors patesiac ec 


een” EOL GRY Wt. (m) 


From the table, the conditions that best fit the values 


used in the numerical calculations are 


Dry density = 250 Keyan? 
Moisture content - 2502 
Depth of peat = Os 25) mt 


The corresponding U-factor values are 


= 1.28 Watt /m 2 
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U-factor in Summer and Winter 


_U-factor in summer, U, = Ulf 


U-factor in winter, U = a 
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TABLE 21. Values of U, and Us for the year 


NS 


Year De U. watts/m--°C 

0.08 = | 0.40 0.16 
O.27 = oe O26 Ons 
O25 == On2o Ona 
0.33 = - 0.20 
0.42 eZo - = 
0.50 1.28 = - 
0.58 ae 2 = - 
0.67 128 - - 
Osi = = 0.68 
0.83 = = 0.37 
O792 - | 0.68 0.26 
1.00 = 0.40 Oras 


The thermal properties of the snow cover used are 250 Hane 
for the density and 0.13 watt/m-°C for the thermal conductivity. These 
properties are dependent on many factors, such as, temperature, moisture 
content, fresh or settled snow and the method used in determining them, 
as such their values varied over a wide range. For example its density 
ranges from 72 to 400 Coe (40, 41) and its thermal conductivity ranges 


from 0.018 to 0.5 watt/m-°C. This wide range in values indicate 
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variability in the thermal properties of snow. The thermal properties 
of peat as those for the snow cover are also dependent on many factors, 
such as, moisture content, porosity, plant structure, temperature and 
drainage. Peat, because of its nature, has a great ability to soak up 
water many times its own weight. The moisture content as such deter- 
mines the thermal properties of peat to a great extent. The moisture 
content of the peat varies considerably, being less in the more 
decomposed types than with more fibrous types, and it also varies with 
season and the corresponding position of the ground water table. For 
pure peats, when saturated, the water content is more than 600 per cent 
of dry weight, generally varying between 750 and 1500 per cent but 
extremes of 50% and 3000% have been noted. However a small percentage 
of inorganic material in the peat will seriously affect the water 
content value, which will drop sharply as the mineral contamination 
increases. The dry unit weight of the peat also varies considerably 
ranging from 40 to 400 aye depending on the type of peat (39). The 
observed thickness of the peat ranges from 30 cm to more than 6 m in 
the discontinuous permafrost ioe It usually does not exceed 1 m in 


the continuous permafrost zone. 
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APPENDIX IV 


Climatic Conditions 
Solar Short Wave Radiation 


The range of values of short wave radiation (Ro) Lotsa 


latitude of 65°N from 1943-1957 (42). 


TABLE 22. Values of R_ for Norman Wells (65°17'N) 
re) 
Months Year Mean Radiation Mean 


Absolute Deviation 


(races) A, PARES) 
January “POTOS 6.94 2 leit 
February On 7 27.78 eee Li, 
March 0.25 94 
April 0.33 533 
May 0.42 a0 
June 0.50 67 
July 0.58 aro 
August 0.67 ~69 
September Oa yA) oo 
October 0.83 Bs) 
November Qe92 78 
December 1.00 - 
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APPENDIX V 


Clausius Clapeyron equation expresses the relationship 


between the temperature T and the saturated vapor pressure e 


ss* 
The equation is 
ae dT i 
Ecates Le wee (i) 
ss (T+273) 
Integrating equation (a) gives 
oe = — X/(T+273). + Y (ii) 


where X and Y are constants. 


TABLE 23. Values of T and e 


ss 


1/ (T+273) Sat. Vapor Press. 


(4) 
0.00341 23:37 
0.00553 UD EBT 
0.00366 6. 
0.00380 2.86 
0.00388 1.91 
0.00395 1 Opes 
0.00403 0.81 
0.00412 Okelew | 
0.00429 0.19 
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X and Y are obtained by obtaining values of See and 
1/(T+273) from Figure (43) and then substituting them into equation 


(1). ‘The: calculated values for X and Y are 


we= 9900) and’ Y'= 23.4 


hence 
a 9234 97 2900/ (1.4273) 
ss 
The air vapor pressure is given by 
e = we 
a ss 
or 


eee h £-3900/ (144273) 


(a0) 
i} 


where w is the relative humidity. 


The surface vapor pressure is given by 
e =M_e + (1-—M )e 
s ss Ss 238 


where M. is the moisture content. 
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FIGURE 43 — Saturated vapor pressure (mb) vs 1/T IC") 
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Month 


January 
February 
March 
April 
May 

June 
July 
August 
September 
October 
November 


December 


TABLE 24. 
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Mean Relative Humidity for 


Norman Wells from Climatic Normals (48) 


Considering only the 
relative humidity w in 
the summer months. The 
average relative humidity 
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APPENDIX VI 


TABLE 25. Mean Monthly Cloud Cover for 
Norman Wells from 


Climatic Normals (49) 


Months 


January 
February 
March 
April 
May 

June 
July 
August 
September 
October 
November 


December 


The mean annual cloudiness We Tomo. 
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TABLE 26. Mean Monthly Wind Velocity (m/sec) 


at Norman Wells from Burns (30) 


Months . ~Mean Monthly Standard Deviation 
Wind Vel. (m/sec) (m/sec) 
January 3.00 Bis20 
February Peet Ideas lye 
March Seon Pls ops 
April Sm Sis we A 
May 3.80 Deol 
June | 4.07 2.41 
July STE: Zio 
August 3.62 2.46 
September : SR UA 2.47 
October S07 2.86 
November 3.08 | 2.02 
December 2003 3.04 


Mean annual wind velocity is 3.43 m/sec. 
Mean annual standard deviation is 2.82 m/sec. 


The value used in the calculations is u = 4.5 m/sec. 
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TABLE 27. Mean Sea-Level Pressure P (mb) for 


Norman Wells from Climatic Normals (49) 


Months P (Mb) 
January LOZ 
February 1024 
March 1022 
April 1020 
May 1018 
June 1014 
July 1013 
August 1013 
September 1014 
October 1012 
November 1018 
pete isi 1a A ee 1018 


ee ood 


Annual mean sea-level pressure P is 1018 mb. This 


corresponds to a ture air pressure at Normal Wells of 1000 mb. 
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